{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b8eefd91-c180-4791-af93-2b24c217257d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from scipy.stats import mannwhitneyu, spearmanr\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests\n",
    "\n",
    "sns.set('talk')\n",
    "sns.set_style('whitegrid')\n",
    "\n",
    "def p_adjust(pvalues, method='fdr_bh'):\n",
    "    res = multipletests(pvalues, method=method)\n",
    "    return np.array(res[1], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "840135a8-d66d-4c69-899b-f283e432947e",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Diversity assocations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f56e9694-1a87-47cb-a219-1539b1d23057",
   "metadata": {},
   "source": [
    "## Read in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1017910f-ba7e-4156-a510-6c2926d4caef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SubmissionType</th>\n",
       "      <th>SampleNumber</th>\n",
       "      <th>SampleIDValidation</th>\n",
       "      <th>DiversigenCheckInSampleName</th>\n",
       "      <th>ReplacesLowVolumeSampleID</th>\n",
       "      <th>BoxLocation</th>\n",
       "      <th>SampleType</th>\n",
       "      <th>SampleSource</th>\n",
       "      <th>SequencingType</th>\n",
       "      <th>BabyN</th>\n",
       "      <th>...</th>\n",
       "      <th>median_mmNorm</th>\n",
       "      <th>median_mmNorm_DTAPHib</th>\n",
       "      <th>median_mmNorm_PCV</th>\n",
       "      <th>PT_protected</th>\n",
       "      <th>Dip_protected</th>\n",
       "      <th>FHA_protected</th>\n",
       "      <th>PRN_protected</th>\n",
       "      <th>TET_protected</th>\n",
       "      <th>PRP (Hib)_protected</th>\n",
       "      <th>VR_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>233_A2_NS_A2_Extra</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1004.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>233_A2_NS_A1</td>\n",
       "      <td>Box 10_D8</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>233.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.231020</td>\n",
       "      <td>0.213158</td>\n",
       "      <td>0.231020</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V3_NS_A2</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1005.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>107_V3_N!_A1</td>\n",
       "      <td>Box 10_D9</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>107.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Not Measured</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>209_V4_NS_A2</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1006.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>209_V4_NS_A1</td>\n",
       "      <td>Box 10_E1</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>209.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.068419</td>\n",
       "      <td>0.035388</td>\n",
       "      <td>0.075616</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V4_NS_2a</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1007.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box9_E2</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>107.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Not Measured</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>229_V2_NS_A1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1008.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 10_C8</td>\n",
       "      <td>Nasal Swab</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>16S</td>\n",
       "      <td>229.0</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Not Measured</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 77 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                     SubmissionType  SampleNumber  SampleIDValidation  \\\n",
       "SampleID                                                                \n",
       "233_A2_NS_A2_Extra  Primary in Tube        1004.0                 NaN   \n",
       "107_V3_NS_A2        Primary in Tube        1005.0                 NaN   \n",
       "209_V4_NS_A2        Primary in Tube        1006.0                 NaN   \n",
       "107_V4_NS_2a        Primary in Tube        1007.0                 NaN   \n",
       "229_V2_NS_A1        Primary in Tube        1008.0                 NaN   \n",
       "\n",
       "                   DiversigenCheckInSampleName ReplacesLowVolumeSampleID  \\\n",
       "SampleID                                                                   \n",
       "233_A2_NS_A2_Extra                         NaN              233_A2_NS_A1   \n",
       "107_V3_NS_A2                               NaN              107_V3_N!_A1   \n",
       "209_V4_NS_A2                               NaN              209_V4_NS_A1   \n",
       "107_V4_NS_2a                               NaN                       NaN   \n",
       "229_V2_NS_A1                               NaN                       NaN   \n",
       "\n",
       "                   BoxLocation  SampleType  SampleSource SequencingType  \\\n",
       "SampleID                                                                  \n",
       "233_A2_NS_A2_Extra   Box 10_D8  Nasal Swab  Human Infant            16S   \n",
       "107_V3_NS_A2         Box 10_D9  Nasal Swab  Human Infant            16S   \n",
       "209_V4_NS_A2         Box 10_E1  Nasal Swab  Human Infant            16S   \n",
       "107_V4_NS_2a           Box9_E2  Nasal Swab  Human Infant            16S   \n",
       "229_V2_NS_A1         Box 10_C8  Nasal Swab  Human Infant            16S   \n",
       "\n",
       "                    BabyN  ...  median_mmNorm median_mmNorm_DTAPHib  \\\n",
       "SampleID                   ...                                        \n",
       "233_A2_NS_A2_Extra  233.0  ...       0.231020              0.213158   \n",
       "107_V3_NS_A2        107.0  ...            NaN                   NaN   \n",
       "209_V4_NS_A2        209.0  ...       0.068419              0.035388   \n",
       "107_V4_NS_2a        107.0  ...            NaN                   NaN   \n",
       "229_V2_NS_A1        229.0  ...            NaN                   NaN   \n",
       "\n",
       "                    median_mmNorm_PCV PT_protected  Dip_protected  \\\n",
       "SampleID                                                            \n",
       "233_A2_NS_A2_Extra           0.231020         True           True   \n",
       "107_V3_NS_A2                      NaN          NaN            NaN   \n",
       "209_V4_NS_A2                 0.075616         True           True   \n",
       "107_V4_NS_2a                      NaN          NaN            NaN   \n",
       "229_V2_NS_A1                      NaN          NaN            NaN   \n",
       "\n",
       "                   FHA_protected  PRN_protected TET_protected  \\\n",
       "SampleID                                                        \n",
       "233_A2_NS_A2_Extra          True           True          True   \n",
       "107_V3_NS_A2                 NaN            NaN           NaN   \n",
       "209_V4_NS_A2               False           True          True   \n",
       "107_V4_NS_2a                 NaN            NaN           NaN   \n",
       "229_V2_NS_A1                 NaN            NaN           NaN   \n",
       "\n",
       "                   PRP (Hib)_protected      VR_group  \n",
       "SampleID                                              \n",
       "233_A2_NS_A2_Extra                True           NVR  \n",
       "107_V3_NS_A2                       NaN  Not Measured  \n",
       "209_V4_NS_A2                      True           NVR  \n",
       "107_V4_NS_2a                       NaN  Not Measured  \n",
       "229_V2_NS_A1                       NaN  Not Measured  \n",
       "\n",
       "[5 rows x 77 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nasal_meta = pd.read_csv('../data/metadata/nasal/nasal_metadata.csv', index_col='SampleID')\n",
    "nasal_meta['age_at_collection'] = (pd.to_datetime(nasal_meta['CollectionDate']) - pd.to_datetime(nasal_meta['DOB'])).dt.days\n",
    "nasal_meta = pd.concat([nasal_meta,\n",
    "                        pd.read_csv('../data/nasal/otu_alpha_diversity.csv', index_col='SampleID'),\n",
    "                        pd.read_csv('../data/nasal/otu_nmds_babies.csv', index_col='SampleID'),\n",
    "                        pd.read_csv('../data/metadata/nasal/nasal_titers_yr2.csv', index_col='SampleID')],\n",
    "                       axis=1)\n",
    "nasal_meta['VR_group'] = nasal_meta['VR_group'].fillna('Not Measured')\n",
    "nasal_meta = nasal_meta.dropna(subset=['VisitCode']) # this drops control samples since they don't have VisitCodes (e.g. MSA2002 and GD5)\n",
    "nasal_meta.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ab9a21c1-9944-49fd-89d4-14627a8bc77b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>SubmissionType</th>\n",
       "      <th>SampleNumber</th>\n",
       "      <th>SampleIDValidation</th>\n",
       "      <th>DiversigenCheckInSampleName</th>\n",
       "      <th>BoxLocation</th>\n",
       "      <th>SampleType</th>\n",
       "      <th>SampleSource</th>\n",
       "      <th>SequencingType</th>\n",
       "      <th>BabyN</th>\n",
       "      <th>BabyN_checked</th>\n",
       "      <th>...</th>\n",
       "      <th>median_mmNorm</th>\n",
       "      <th>median_mmNorm_DTAPHib</th>\n",
       "      <th>median_mmNorm_PCV</th>\n",
       "      <th>PT_protected</th>\n",
       "      <th>Dip_protected</th>\n",
       "      <th>FHA_protected</th>\n",
       "      <th>PRN_protected</th>\n",
       "      <th>TET_protected</th>\n",
       "      <th>PRP (Hib)_protected</th>\n",
       "      <th>VR_group</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SampleID</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>204_V5</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>204_S</td>\n",
       "      <td>Box 7, A1</td>\n",
       "      <td>Stool</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>MetaG</td>\n",
       "      <td>204.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.035342</td>\n",
       "      <td>0.027508</td>\n",
       "      <td>0.060708</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>226_V1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 7, A2</td>\n",
       "      <td>Stool</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>MetaG</td>\n",
       "      <td>226.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.082202</td>\n",
       "      <td>0.097174</td>\n",
       "      <td>0.057712</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107_V3</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>3.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 7, A3</td>\n",
       "      <td>Stool</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>MetaG</td>\n",
       "      <td>107.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>108_V3</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>4.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 7, A4</td>\n",
       "      <td>Stool</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>MetaG</td>\n",
       "      <td>108.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>0.168430</td>\n",
       "      <td>0.172291</td>\n",
       "      <td>0.164958</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>NVR</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>109_V1</th>\n",
       "      <td>Primary in Tube</td>\n",
       "      <td>5.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>Box 7, A5</td>\n",
       "      <td>Stool</td>\n",
       "      <td>Human Infant</td>\n",
       "      <td>MetaG</td>\n",
       "      <td>109.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 88 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "           SubmissionType  SampleNumber  SampleIDValidation  \\\n",
       "SampleID                                                      \n",
       "204_V5    Primary in Tube           1.0                 NaN   \n",
       "226_V1    Primary in Tube           2.0                 NaN   \n",
       "107_V3    Primary in Tube           3.0                 NaN   \n",
       "108_V3    Primary in Tube           4.0                 NaN   \n",
       "109_V1    Primary in Tube           5.0                 NaN   \n",
       "\n",
       "         DiversigenCheckInSampleName BoxLocation SampleType  SampleSource  \\\n",
       "SampleID                                                                    \n",
       "204_V5                         204_S   Box 7, A1      Stool  Human Infant   \n",
       "226_V1                           NaN   Box 7, A2      Stool  Human Infant   \n",
       "107_V3                           NaN   Box 7, A3      Stool  Human Infant   \n",
       "108_V3                           NaN   Box 7, A4      Stool  Human Infant   \n",
       "109_V1                           NaN   Box 7, A5      Stool  Human Infant   \n",
       "\n",
       "         SequencingType  BabyN  BabyN_checked  ... median_mmNorm  \\\n",
       "SampleID                                       ...                 \n",
       "204_V5            MetaG  204.0            NaN  ...      0.035342   \n",
       "226_V1            MetaG  226.0            NaN  ...      0.082202   \n",
       "107_V3            MetaG  107.0            NaN  ...           NaN   \n",
       "108_V3            MetaG  108.0            NaN  ...      0.168430   \n",
       "109_V1            MetaG  109.0            NaN  ...           NaN   \n",
       "\n",
       "          median_mmNorm_DTAPHib median_mmNorm_PCV  PT_protected Dip_protected  \\\n",
       "SampleID                                                                        \n",
       "204_V5                 0.027508          0.060708         False          True   \n",
       "226_V1                 0.097174          0.057712          True          True   \n",
       "107_V3                      NaN               NaN           NaN           NaN   \n",
       "108_V3                 0.172291          0.164958          True          True   \n",
       "109_V1                      NaN               NaN           NaN           NaN   \n",
       "\n",
       "         FHA_protected  PRN_protected TET_protected  PRP (Hib)_protected  \\\n",
       "SampleID                                                                   \n",
       "204_V5           False           True          True                 True   \n",
       "226_V1            True           True          True                 True   \n",
       "107_V3             NaN            NaN           NaN                  NaN   \n",
       "108_V3            True           True          True                 True   \n",
       "109_V1             NaN            NaN           NaN                  NaN   \n",
       "\n",
       "          VR_group  \n",
       "SampleID            \n",
       "204_V5         NVR  \n",
       "226_V1         NVR  \n",
       "107_V3         NaN  \n",
       "108_V3         NVR  \n",
       "109_V1         NaN  \n",
       "\n",
       "[5 rows x 88 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stool_nmds = pd.read_csv('../data/stool/kraken_nmds_babies.csv', index_col='SampleID')\n",
    "stool_nmds.columns = ['kraken_' + i if i.startswith('MDS') else i for i in stool_nmds.columns]\n",
    "stool_ko_nmds = pd.read_csv('../data/stool/ko_nmds_babies.csv', index_col='SampleID')\n",
    "stool_ko_nmds.columns = ['ko_' + i if i.startswith('MDS') else i for i in stool_ko_nmds.columns]\n",
    "stool_meta = pd.concat([pd.read_csv('../data/metadata/stool/stool_metadata.csv', index_col='SampleID'),\n",
    "                        pd.read_csv('../data/stool/kraken_alpha_diversity.csv', index_col='SampleID'),\n",
    "                        stool_nmds,\n",
    "                        pd.read_csv('../data/stool/ko_alpha_diversity.csv', index_col='SampleID'),\n",
    "                        stool_ko_nmds,\n",
    "                        pd.read_csv('../data/metadata/stool/stool_titers_yr2.csv', index_col='SampleID')],\n",
    "                       axis=1)\n",
    "stool_meta = stool_meta.dropna(subset=['VisitCode'])\n",
    "stool_meta = stool_meta.query('`gt_2.5`')\n",
    "stool_meta.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3458549-3c60-4301-a41c-fb1a98b3bf3d",
   "metadata": {},
   "source": [
    "## Nasal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5e2c92a5-5c11-441d-a954-adcdd72af59d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>stat</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [timepoint, div_metric, titer_measure, N_LVR, N_NVR, stat, p_value]\n",
       "Index: []"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nasal_titer_assoc_rows = list()\n",
    "# for timepoint in ['V5', 'V6', 'V7', 'V9']:\n",
    "for timepoint in set([i for i in nasal_meta['VisitCode'] if i.startswith('V')]):\n",
    "    for div_metric in ['n_otus_div', 'simpson_e_div', 'shannon_div', 'MDS1', 'MDS2']:\n",
    "        N_LVR = len(nasal_meta.query(\"VisitCode == @timepoint & VR_group == 'LVR'\")[div_metric].dropna())\n",
    "        N_NVR = len(nasal_meta.query(\"VisitCode == @timepoint & VR_group == 'NVR'\")[div_metric].dropna())\n",
    "        if N_LVR > 2 and N_NVR > 2:\n",
    "            w_value, p_value = mannwhitneyu(nasal_meta.query(\"VisitCode == @timepoint & VR_group == 'LVR'\")[div_metric].values,\n",
    "                                            nasal_meta.query(\"VisitCode == @timepoint & VR_group == 'NVR'\")[div_metric].values,\n",
    "                                           nan_policy='omit')\n",
    "            nasal_titer_assoc_rows.append([timepoint, div_metric, 'LVR/NVR', N_LVR, N_NVR, w_value, p_value])\n",
    "nasal_titer_assocs = pd.DataFrame(nasal_titer_assoc_rows, columns=['timepoint', 'div_metric', 'titer_measure', 'N_LVR', 'N_NVR', 'stat', 'p_value']).sort_values('p_value')\n",
    "nasal_titer_assocs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c84a0c87-c61a-4fda-836b-71bb9bb762a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>50</td>\n",
       "      <td>0.454801</td>\n",
       "      <td>0.000906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>55</td>\n",
       "      <td>-0.387002</td>\n",
       "      <td>0.003841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>50</td>\n",
       "      <td>0.329268</td>\n",
       "      <td>0.019548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>55</td>\n",
       "      <td>-0.313894</td>\n",
       "      <td>0.020810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td>V6</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>52</td>\n",
       "      <td>0.293825</td>\n",
       "      <td>0.045005</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    timepoint     div_metric          titer_measure   N         R   p_value\n",
       "145        V8           MDS1  median_mmNorm_DTAPHib  50  0.454801  0.000906\n",
       "95         V9  simpson_e_div      median_mmNorm_PCV  55 -0.387002  0.003841\n",
       "144        V8           MDS1          median_mmNorm  50  0.329268  0.019548\n",
       "93         V9  simpson_e_div          median_mmNorm  55 -0.313894  0.020810\n",
       "129        V6           MDS1          median_mmNorm  52  0.293825  0.045005"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nasal_titer_correls_rows = list()\n",
    "# for timepoint in ['V5', 'V6', 'V7', 'V9']:\n",
    "for timepoint in set([i for i in nasal_meta['VisitCode'] if i.startswith('V')]):\n",
    "    for div_metric in ['n_otus_div', 'simpson_e_div', 'shannon_div', 'MDS1', 'MDS2']:\n",
    "        for titer_measure in ['median_mmNorm', 'median_mmNorm_DTAPHib', 'median_mmNorm_PCV']:\n",
    "            N = len(nasal_meta.query(\"VisitCode == @timepoint\")[titer_measure].dropna())\n",
    "            if N > 3:\n",
    "                r_value, p_value = spearmanr(nasal_meta.query(\"VisitCode == @timepoint\")[div_metric].values,\n",
    "                                             nasal_meta.query(\"VisitCode == @timepoint\")[titer_measure].values,\n",
    "                                             nan_policy='omit')\n",
    "                nasal_titer_correls_rows.append([timepoint, div_metric, titer_measure, N, r_value, p_value])\n",
    "nasal_titer_correls = pd.DataFrame(nasal_titer_correls_rows, columns=['timepoint', 'div_metric', 'titer_measure', 'N', 'R', 'p_value']).sort_values('p_value')\n",
    "nasal_titer_correls.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6e7acb48-9ab5-4c73-a0ca-f948e0f87492",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.454801</td>\n",
       "      <td>0.000906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.387002</td>\n",
       "      <td>0.003841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.329268</td>\n",
       "      <td>0.019548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.313894</td>\n",
       "      <td>0.02081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V6</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>52.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.293825</td>\n",
       "      <td>0.045005</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint     div_metric          titer_measure N_LVR N_NVR  \\\n",
       "145     nasal        V8           MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "95      nasal        V9  simpson_e_div      median_mmNorm_PCV   NaN   NaN   \n",
       "144     nasal        V8           MDS1          median_mmNorm   NaN   NaN   \n",
       "93      nasal        V9  simpson_e_div          median_mmNorm   NaN   NaN   \n",
       "129     nasal        V6           MDS1          median_mmNorm   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value  \n",
       "145  50.0  NaN  0.454801  0.000906  \n",
       "95   55.0  NaN -0.387002  0.003841  \n",
       "144  50.0  NaN  0.329268  0.019548  \n",
       "93   55.0  NaN -0.313894   0.02081  \n",
       "129  52.0  NaN  0.293825  0.045005  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nasal_col_order = ['body_site', 'timepoint', 'div_metric', 'titer_measure', 'N_LVR', 'N_NVR', 'N', 'stat', 'R', 'p_value']\n",
    "nasal_titer_stats = pd.concat([nasal_titer_assocs, nasal_titer_correls]).sort_values('p_value')\n",
    "nasal_titer_stats['body_site'] = 'nasal'\n",
    "nasal_titer_stats = nasal_titer_stats[nasal_col_order]\n",
    "nasal_titer_stats.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e930c32c-d338-47c2-8bc1-d589bd119c25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.454801</td>\n",
       "      <td>0.000906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.387002</td>\n",
       "      <td>0.003841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.329268</td>\n",
       "      <td>0.019548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.313894</td>\n",
       "      <td>0.02081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V6</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>52.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.293825</td>\n",
       "      <td>0.045005</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint     div_metric          titer_measure N_LVR N_NVR  \\\n",
       "145     nasal        V8           MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "95      nasal        V9  simpson_e_div      median_mmNorm_PCV   NaN   NaN   \n",
       "144     nasal        V8           MDS1          median_mmNorm   NaN   NaN   \n",
       "93      nasal        V9  simpson_e_div          median_mmNorm   NaN   NaN   \n",
       "129     nasal        V6           MDS1          median_mmNorm   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value  \n",
       "145  50.0  NaN  0.454801  0.000906  \n",
       "95   55.0  NaN -0.387002  0.003841  \n",
       "144  50.0  NaN  0.329268  0.019548  \n",
       "93   55.0  NaN -0.313894   0.02081  \n",
       "129  52.0  NaN  0.293825  0.045005  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nasal_titer_stats.query('p_value < .05')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a52b641-c734-464d-b867-d495db1ca82d",
   "metadata": {},
   "source": [
    "## Stool"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2f7b8c7d-5624-4c5e-b3ec-ea9512e53060",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>stat</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [timepoint, div_metric, titer_measure, N_LVR, N_NVR, stat, p_value]\n",
       "Index: []"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stool_titer_assoc_rows = list()\n",
    "# for timepoint in ['V5', 'V6', 'V7', 'V9']:\n",
    "for timepoint in set([i for i in stool_meta['VisitCode'] if i.startswith('V')]):\n",
    "    for div_metric in ['genus_richness', 'genus_evenness', 'genus_shannon', 'ko_richness', 'ko_evenness', 'ko_shannon', 'kraken_MDS1', 'kraken_MDS2', 'ko_MDS1', 'ko_MDS2']:\n",
    "        N_LVR = len(stool_meta.query(\"VisitCode == @timepoint & VR_group == 'LVR'\")[div_metric].dropna())\n",
    "        N_NVR = len(stool_meta.query(\"VisitCode == @timepoint & VR_group == 'NVR'\")[div_metric].dropna())\n",
    "        if N_LVR > 2 and N_NVR > 2:\n",
    "            w_value, p_value = mannwhitneyu(stool_meta.query(\"VisitCode == @timepoint & VR_group == 'LVR'\")[div_metric].values,\n",
    "                                            stool_meta.query(\"VisitCode == @timepoint & VR_group == 'NVR'\")[div_metric].values,\n",
    "                                            nan_policy='omit')\n",
    "            stool_titer_assoc_rows.append([timepoint, div_metric, 'LVR/NVR', N_LVR, N_NVR, w_value, p_value])\n",
    "stool_titer_assocs = pd.DataFrame(stool_titer_assoc_rows, columns=['timepoint', 'div_metric', 'titer_measure', 'N_LVR', 'N_NVR', 'stat', 'p_value']).sort_values('p_value')\n",
    "stool_titer_assocs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1198d99d-1ed8-4844-8489-b88bf3201148",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>341</th>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>20</td>\n",
       "      <td>-0.556600</td>\n",
       "      <td>0.010806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>303</th>\n",
       "      <td>V2</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>26</td>\n",
       "      <td>-0.478974</td>\n",
       "      <td>0.013302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>339</th>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>20</td>\n",
       "      <td>-0.522001</td>\n",
       "      <td>0.018233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>328</th>\n",
       "      <td>V2</td>\n",
       "      <td>ko_MDS2</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>26</td>\n",
       "      <td>-0.459145</td>\n",
       "      <td>0.018299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>V4</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>14</td>\n",
       "      <td>0.604396</td>\n",
       "      <td>0.022057</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    timepoint      div_metric          titer_measure   N         R   p_value\n",
       "341        V3     ko_richness      median_mmNorm_PCV  20 -0.556600  0.010806\n",
       "303        V2  genus_evenness          median_mmNorm  26 -0.478974  0.013302\n",
       "339        V3     ko_richness          median_mmNorm  20 -0.522001  0.018233\n",
       "328        V2         ko_MDS2  median_mmNorm_DTAPHib  26 -0.459145  0.018299\n",
       "150        V4  genus_richness          median_mmNorm  14  0.604396  0.022057"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stool_titer_correls_rows = list()\n",
    "# for timepoint in ['V5', 'V6', 'V7', 'V9']:\n",
    "for timepoint in set([i for i in stool_meta['VisitCode'] if i.startswith('V')]):\n",
    "    for div_metric in ['genus_richness', 'genus_evenness', 'genus_shannon', 'ko_richness', 'ko_evenness', 'ko_shannon', 'kraken_MDS1', 'kraken_MDS2', 'ko_MDS1', 'ko_MDS2']:\n",
    "        for titer_measure in ['median_mmNorm', 'median_mmNorm_DTAPHib', 'median_mmNorm_PCV']:\n",
    "            N = len(stool_meta.query(\"VisitCode == @timepoint\")[titer_measure].dropna())\n",
    "            if N > 3:\n",
    "                r_value, p_value = spearmanr(stool_meta.query(\"VisitCode == @timepoint\")[div_metric].values,\n",
    "                                             stool_meta.query(\"VisitCode == @timepoint\")[titer_measure].values,\n",
    "                                             nan_policy='omit')\n",
    "                stool_titer_correls_rows.append([timepoint, div_metric, titer_measure, N, r_value, p_value])\n",
    "stool_titer_correls = pd.DataFrame(stool_titer_correls_rows, columns=['timepoint', 'div_metric', 'titer_measure', 'N', 'R', 'p_value']).sort_values('p_value')\n",
    "stool_titer_correls.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "242e75e1-c1bf-4b62-aedb-db64dfc3137f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>341</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.556600</td>\n",
       "      <td>0.010806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>303</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.478974</td>\n",
       "      <td>0.013302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>339</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.522001</td>\n",
       "      <td>0.018233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>328</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>ko_MDS2</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.459145</td>\n",
       "      <td>0.018299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>stool</td>\n",
       "      <td>V4</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.604396</td>\n",
       "      <td>0.022057</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint      div_metric          titer_measure N_LVR N_NVR  \\\n",
       "341     stool        V3     ko_richness      median_mmNorm_PCV   NaN   NaN   \n",
       "303     stool        V2  genus_evenness          median_mmNorm   NaN   NaN   \n",
       "339     stool        V3     ko_richness          median_mmNorm   NaN   NaN   \n",
       "328     stool        V2         ko_MDS2  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "150     stool        V4  genus_richness          median_mmNorm   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value  \n",
       "341  20.0  NaN -0.556600  0.010806  \n",
       "303  26.0  NaN -0.478974  0.013302  \n",
       "339  20.0  NaN -0.522001  0.018233  \n",
       "328  26.0  NaN -0.459145  0.018299  \n",
       "150  14.0  NaN  0.604396  0.022057  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stool_col_order = ['body_site', 'timepoint', 'div_metric', 'titer_measure', 'N_LVR', 'N_NVR', 'N', 'stat', 'R', 'p_value']\n",
    "stool_titer_stats = pd.concat([stool_titer_assocs, stool_titer_correls]).sort_values('p_value')\n",
    "stool_titer_stats['body_site'] = 'stool'\n",
    "stool_titer_stats = stool_titer_stats[stool_col_order]\n",
    "stool_titer_stats.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "548f02c2-f678-4102-afba-af83e3bdaa6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>341</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.556600</td>\n",
       "      <td>0.010806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>303</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.478974</td>\n",
       "      <td>0.013302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>339</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.522001</td>\n",
       "      <td>0.018233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>328</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>ko_MDS2</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.459145</td>\n",
       "      <td>0.018299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>stool</td>\n",
       "      <td>V4</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.604396</td>\n",
       "      <td>0.022057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.352104</td>\n",
       "      <td>0.032579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>349</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.470677</td>\n",
       "      <td>0.036215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.342579</td>\n",
       "      <td>0.037933</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.450649</td>\n",
       "      <td>0.040345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>genus_shannon</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.448052</td>\n",
       "      <td>0.041656</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>stool</td>\n",
       "      <td>V5</td>\n",
       "      <td>ko_evenness</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>32.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.354106</td>\n",
       "      <td>0.046769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.433766</td>\n",
       "      <td>0.049467</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint      div_metric          titer_measure N_LVR N_NVR  \\\n",
       "341     stool        V3     ko_richness      median_mmNorm_PCV   NaN   NaN   \n",
       "303     stool        V2  genus_evenness          median_mmNorm   NaN   NaN   \n",
       "339     stool        V3     ko_richness          median_mmNorm   NaN   NaN   \n",
       "328     stool        V2         ko_MDS2  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "150     stool        V4  genus_richness          median_mmNorm   NaN   NaN   \n",
       "180     stool        V9  genus_richness          median_mmNorm   NaN   NaN   \n",
       "349     stool        V3     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "199     stool        V9     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "95      stool        V1  genus_evenness      median_mmNorm_PCV   NaN   NaN   \n",
       "97      stool        V1   genus_shannon  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "43      stool        V5     ko_evenness  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "110     stool        V1     kraken_MDS1      median_mmNorm_PCV   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value  \n",
       "341  20.0  NaN -0.556600  0.010806  \n",
       "303  26.0  NaN -0.478974  0.013302  \n",
       "339  20.0  NaN -0.522001  0.018233  \n",
       "328  26.0  NaN -0.459145  0.018299  \n",
       "150  14.0  NaN  0.604396  0.022057  \n",
       "180  37.0  NaN  0.352104  0.032579  \n",
       "349  20.0  NaN  0.470677  0.036215  \n",
       "199  37.0  NaN  0.342579  0.037933  \n",
       "95   21.0  NaN -0.450649  0.040345  \n",
       "97   21.0  NaN  0.448052  0.041656  \n",
       "43   32.0  NaN -0.354106  0.046769  \n",
       "110  21.0  NaN  0.433766  0.049467  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stool_titer_stats.query('p_value < .05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4aa67a50-5300-4628-948d-4741af77cc87",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Meidan DTAPHib Titer')"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHHCAYAAACfqw0dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB0nUlEQVR4nO3deVxU1f8/8NcwDDCA4oKAgDsOCuLHHUFzjUwzcyFb3DJLs800s0XbLLePaWZq2mKFW6loiZr4SdCv5IIbuSUKoiIqEAKKw8Aw3N8f/GaCGODOZQZm4PV8PHiU95x77pnDZeY955x7jkwQBAFEREREZBK72q4AERERkS1iEEVEREQkAYMoIiIiIgkYRBERERFJwCCKiIiISAIGUUREREQSMIgiIiIiksC+titQl507dw5arRZ2dnZwdHSs7eoQERGRCAUFBSguLoZCoUBQUFCF+RhEWZBWq4UgCNDpdFCr1bVdHSIiIjKBVqutNJ1BlAXZ2dlBp9NBJpNBqVRWmE8QBOTn5wMAlEolZDJZTVWx3mAbWxbb17LYvpbF9rU8W2vj/Px8CIIAO7vKZz0xiLIgR0dHqNVqKJVKdOzYscJ8Op0OCQkJAAB/f3/I5fIaqmH9wTa2LLavZbF9LYvta3m21sZ//fUX1Gp1lVNxOLGciIiISAKb6YlKSUnB6tWrcerUKWRlZcHLywtDhw7FtGnT4OzsLKqMCRMmID4+vsp8Pj4+iImJqW6ViYiIqA6ziSDq7NmzmDRpEtRqNTp37oygoCCcPn0aa9euRWxsLDZv3gxXV9cqywkNDYWnp2eF6bGxscjLy0NgYKA5q09ERER1kNUHUUVFRZg1axbUajUWLFiA8PBwAIBGo8HMmTMRExOD5cuX44MPPqiyrOnTp1eY9ssvvyAqKgpt2rTBokWLzFZ/IiIiqpusfk7Unj17kJqaipCQEEMABQBOTk5YuHAhnJ2dsXXrVuTm5kq+RnJyMj766CMoFAp8/vnnonq1iIiIqH6z+iBKPzcpLCysXFrjxo0RHBwMrVaLw4cPS77Gxx9/jPz8fEyZMqXSp+iIiIiI9Kx+OO/y5csASh6JNMbPzw+xsbG4dOkShg8fbnL5e/bswfHjx9G8eXO89NJL1aprRfQLblakdFpl+Ug6trFlsX0ti+1rWWxfy7O1NhYEQVQ+qw+iMjIyAKDCCeEeHh5l8pmiuLgYX3zxBQBg2rRplS6IWR35+fmG9TGqcu7cOYvUwZYolUrI5XLodDrD4mzmxDa2LLavZbF9LYvta3l1qY2tPojSb5fi5ORkNF1/XMq2KtHR0bh+/To8PDwwZswY6ZWkapPL5WjatBkaN20KtUaH++pCNHZ2gLOjHNl3s5CVlWkT316IiKj+sPogSi6Xo7i4uMol4sV2vZX2ww8/AACee+45ODg4SKmeKEqlssLhSKCka1MfmQcFBVn9Sq6WoC0SEH/xDj7bdgxX0/55SKCtjxvGDPRDr4BOUNhL3yaAbWxZbF/LYvtaFtvX8mytjRMTE0WNhFh9EOXi4oKcnJwKX4xGowEAk4fiUlNTkZCQALlcjieeeKLa9ayMTCYTfcPI5XKrv7nMTVNYhPW7L+C3I9fKpV1Ny8XSjacwLLQ1Jj8eCCeH6t+y9bGNaxLb17LYvpbF9rU8W2hjsXv7Wf3Tefo5T5mZmUbT9XOh9PnE2rdvHwCgd+/ecHd3r0YNqTp0umLEX7hjNIAqbe+RazhxIR06XXHNVIyIiKgKVh9E6YfBkpKSjKbrj1c2XGbMoUOHAACPPvpoNWpH1aXVFWN7zBVReSNjr0DLIIqIiKyE1QdR/fv3B1AyCfzfsrOzcfz4cSgUCvTp00d0mTqdDufPnwcA9OjRwzwVJUke5GuRcuueqLzJabl4kK+1cI2IiIjEsfogKiwsDN7e3oiLi8OmTZsMxzUaDebOnQu1Wo3w8PAyQ3JarRbJyclITk6GVlv+Q/fKlSvIz89Ho0aN0LZt2xp5HWRcntq0oCiPQRQREVkJq59Y7uTkhMWLF2Pq1KmYP38+IiMj4evrizNnziAjIwMBAQGYPXt2mXPS09MxbNgwAMCBAwfg6+tbJv3mzZsAUO441TxXZ4Vp+ZWm5SciIrIUq++JAoDg4GBs27YNQ4YMwa1bt3Dw4EE0aNAAL7/8MjZs2GDyXnd3794FADRv3twS1SUTuCgVaOvjJipvOx83uDgxiCIiIutg9T1ReiqVCitXrhSV19fXF4mJiRWmjx07FmPHjjVX1agaFHI7jBnoh6UbT1WZd8zA9lDY20TcT0RE9QA/kahWyeV26BXohaGhrSvNNyy0NXoGekIu5y1LRETWwWZ6oqjucnKwx/OPByKorTu2x14pt2J5+MD26BnoaZaFNomIiMyFn0pkFZwc7BHauTl6BnriQb4WeflauCoVcHFSQGFvxx4oIiKyOgyiyGrI5SXBkpODPZq6mbaNDxERUU3j13siIiIiCRhEEREREUnAIIqIiIhIAgZRRERERBIwiCIiIiKSgEEUERERkQQMooiIiIgkYBBFREREJAGDKCIiIiIJGEQRERERScAgioiIiEgCBlFEREREEjCIIiIiIpKAQRQRERGRBAyiiIiIiCRgEEVEREQkAYMoIiIiIgkYRBERERFJwCCKiIiISAIGUUREREQSMIgiIiIikoBBFBEREZEEDKKIiIiIJGAQRURERCQBgygiIiIiCRhEEREREUnAIIqIiIhIAgZRRERERBLY13YFxEpJScHq1atx6tQpZGVlwcvLC0OHDsW0adPg7OxsUllqtRrr169HdHQ0bty4ATs7O3Ts2BETJ07Eo48+aqFXQERERHWJTfREnT17FqNHj0ZUVBTc3d0xYMAAqNVqrF27Fk8//TTy8vJEl/X333/jySefxJdffom7d++ib9++8Pf3x+nTpzFjxgz88MMPlnshREREVGdYfRBVVFSEWbNmQa1WY8GCBdi2bRtWrlyJ33//HYMGDUJiYiKWL18uury5c+ciKSkJQ4YMQUxMDFavXo2ffvoJ69evh0KhwH//+1/cunXLgq+IiIiI6gKrD6L27NmD1NRUhISEIDw83HDcyckJCxcuhLOzM7Zu3Yrc3Nwqyzp79iwOHjyIli1bYunSpXB0dDSkhYaGYvTo0fDy8sL58+ct8lqIiIio7rD6IComJgYAEBYWVi6tcePGCA4OhlarxeHDh6ss67fffgMATJgwoUwApTd//nzExMTgkUceqWatiYiIqK6z+onlly9fBgD4+/sbTffz80NsbCwuXbqE4cOHV1qWvoepS5cuUKvViI6Oxrlz56DT6RAUFIThw4fDycnJvC+AiIiI6iSrD6IyMjIAAJ6enkbTPTw8yuSrzLVr1wAA2dnZGD58ONLS0gxpP/30E7766iusXbsW7du3r2atyxIEATqdrsL00mmV5SPp2MaWxfa1LLavZbF9Lc/W2lgQBFH5rD6IUqvVAFBhD5H+uD5fZfRP8c2ePRve3t7YuHEjOnbsiJs3b+Kzzz7D4cOH8eKLL2L37t1wdXU10ysA8vPzkZCQICrvuXPnzHZdMo5tbFlsX8ti+1oW29fy6lIbW/2cKLlcDgCQyWSV5hMTNRYUFAAAHBwcEBERgZ49e8LV1RUdOnTA2rVroVKpcPv2bfz000/VrzgRERHVaZJ7olJSUtC6desqg5vqcnFxQU5ODvLz842mazQaAIBSqayyLCcnJzx48AAjR46Em5tbmTR7e3s8/fTTmD9/Po4ePYoXXnih+pX//5RKZYVzuoCSrk19ZB4UFGQIHMl82MaWxfa1LLavZbF9Lc/W2jgxMbHCuKM0yUHU66+/jgcPHiAyMhKNGzeWWkyVPDw8kJOTg8zMTLRo0aJcun4ulH5uVGXc3d3x4MED+Pr6Gk3XH8/Ozq5GjcuTyWSibxi5XG71N5etYxtbFtvXsti+lsX2tTxbaGOxHUSSh/NSU1Ph6Oho0QAK+OepvKSkJKPp+uOV9fT8u6z09HSj6ZmZmQCAJk2amFxPIiIiql8kB1EuLi4oLi42Z12M6t+/PwAgOjq6XFp2djaOHz8OhUKBPn36VFnWgAEDAJSsF6XVasul/9///R8AoFevXtWoMREREdUHkoOoZ599FtevX8f69evNWZ9ywsLC4O3tjbi4OGzatMlwXKPRYO7cuVCr1QgPD4e7u7shTavVIjk5GcnJyWWCpWHDhsHX1xfXrl3D/Pnzy6Rt27YN0dHRcHNzw5gxYyz6moiIiMj2SZ4TFRQUhG7dumHp0qXYtGkTunbtimbNmlW6WOWMGTNMvo6TkxMWL16MqVOnYv78+YiMjISvry/OnDmDjIwMBAQEYPbs2WXOSU9Px7BhwwAABw4cMMx1UiqV+OKLL/DCCy9g69atOHjwIP7zn//g+vXruHz5MhwdHbF48WI0bdrU5HoSERFR/SI5iJo6dSpkMhkEQUBaWlqlm/YKggCZTCYpiAKA4OBgbNu2DatWrUJ8fDySkpLg6+uL8PBwTJkyxaQ1nTp16oSoqCisW7cOsbGxOHjwIBo1aoTHHnsMU6dORYcOHSTVkYiIiOoXyUFUz549zVmPKqlUKqxcuVJUXl9fXyQmJlaY3qxZM8ybNw/z5s0zV/WIiIionpEcRG3YsMGc9SAiIiKyKVa/YjkRERGRNTLL3nm5ubk4evQorl69ivv37+Ptt99GQUEB/vzzTy4XQERERHVStYIoQRDw5Zdf4vvvvzdsvwIAb7/9Nm7evIlJkyahS5cuWL16NRewJCIiojqlWsN5c+bMwVdffYX8/Hw0atSozP51OTk5EAQBCQkJmDBhgqg9aIiIiIhsheQgav/+/YiKikKTJk3wzTff4OjRo2WWB+jevTs2bdqEpk2b4urVq4iIiDBLhYmIiIisgeQg6ueff4ZMJsOyZcvw0EMPGc3TvXt3rFixAoIgGN22hYiIiMhWSQ6izp8/j+bNm6N3796V5uvRowd8fHxw7do1qZciIiIisjqSgyi1Wo1GjRqJytukSRMUFRVJvRQRERGR1ZEcRLm7u+P69esQBKHSfFqtFteuXSuzQTARERGRrZMcRPXq1QtqtRqbN2+uNN+PP/6I+/fvo0ePHlIvRURERGR1JAdRkydPhp2dHZYsWYKIiAhkZ2eXSc/KysKKFSuwfPly2NnZYfz48dWuLBEREZG1kLzYZocOHfDee+/h008/xaJFi7Bo0SJDWkhICHJycgCULMg5Y8YMdO7cudqVJSIiIrIW1Vpsc9y4cVi3bh38/f0hCILhJzs7G4IgoGXLlli+fDmmT59urvoSERERWYVq753Xr18/9OvXD2lpabhy5Qru378PpVKJNm3aoF27duaoIxEREZHVkRxEnThxAg0aNDCsUu7j4wMfHx+jeePi4nD9+nWMGzdO6uWIiIiIrIrk4bwJEybg008/FZV3xYoVWLFihdRLEREREVkdUT1ReXl55Z6+AwCNRoPU1NQKzxMEAWlpabh69WqV60kRERER2RJRQdSDBw/w+OOPo6CgwHBMJpPhwoULeOSRR0RdqEuXLpIqSERERGSNRA3neXp6YvLkyWWewANQ5t+V/TRv3hzz5s2z6AshIiIiqkmiJ5a//PLLCA8PB1ASPD388MMICgqqdK6TnZ0dnJ2d4ebmVu2KEhEREVkT0UGUQqEo8/Rdz5494e/vX+ETeURERER1meQlDjZs2GDOehARERHZFFFBVHFxMYCS4bl/HzNF6fOJiIiIbJmoICogIAB2dnbYs2cP2rRpAwAIDAw06UIymQwXL140vYZEREREVkj0cN6/e5647hMRERHVZ6KCqIiICACAt7d3uWNERERE9ZGoICo+Ph7e3t7o1auX4Vjp/yciIiKqb0TN9F61ahUiIyMtXRciIiIim8HH5YiIiIgkYBBFREREJAGDKCIiIiIJGEQRERERSSB6najr16/j3XfflXwhmUyGhQsXSj4/JSUFq1evxqlTp5CVlQUvLy8MHToU06ZNg7Ozs+hyNBoNunXrBp1OV2Ge7du3IygoSHJdiYiIqO4THURlZWXhl19+kXQRQRCqFUSdPXsWkyZNglqtRufOnREUFITTp09j7dq1iI2NxebNm+Hq6iqqrEuXLkGn08Hb2xvdu3c3mqdRo0aS6klERET1h+ggytXVFR06dLBkXYwqKirCrFmzoFarsWDBAoSHhwMo6VGaOXMmYmJisHz5cnzwwQeiyrtw4QIA4PHHH8esWbMsVm8iIiKq20QHUSqVChs2bLBkXYzas2cPUlNTERISYgigAMDJyQkLFy7EoEGDsHXrVsyYMQNubm5VlqcPojhcR0RERNVh9RPLY2JiAABhYWHl0ho3bozg4GBotVocPnxYVHkMooiIiMgcrD6Iunz5MgDA39/faLqfnx+AkrlOVSksLERycjIaNWqEo0eP4umnn0aPHj3QvXt3TJ48GX/88Yf5Kk5ERER1mujhvNqSkZEBAPD09DSa7uHhUSZfZS5dugStVoucnBy8++676NKlC4KDg5GUlIQjR47gyJEjePPNNzF16lTzvQCUTKyv7GnA0mmV5SPp2MaWxfa1LLavZbF9Lc/W2lgQBFH5rD6IUqvVAErmQBmjP67PV5mLFy8CKAm81qxZU2ZIb+fOnZg7dy6WL1+Orl27omfPntWtukF+fj4SEhJE5T137pzZrkvGsY0ti+1rWWxfy2L7Wl5damNRQdSoUaPQunVrC1fFOLlcjuLiYshkskrziYkax44di379+kEul5fr2Ro1ahQuXLiADRs2ICIiwqxBFBEREdU9ooKoRYsWWboeFXJxcUFOTg7y8/ONpms0GgCAUqmssiw7Ozt4e3tXmD548GBs2LDB7FGyUqmscE4XUNK1qb9mUFAQ5HK5Wa9PbGNLY/taFtvXsti+lmdrbZyYmFhh3FGa1Q/neXh4ICcnB5mZmWjRokW5dP1cKP3cqOrw8vICAFENZwqZTCb6hpHL5VZ/c9k6trFlsX0ti+1rWWxfy7OFNq5q9EvP6p/O0/fgJCUlGU3XH6+sp0dvzZo1eP3113H06FGj6Xfu3AHwTzBFREREVBGrD6L69+8PAIiOji6Xlp2djePHj0OhUKBPnz5VlpWSkoLo6Gjs3LnTaLr++IABA6RXmIiIiOoFqw+iwsLC4O3tjbi4OGzatMlwXKPRYO7cuVCr1QgPD4e7u7shTavVIjk5GcnJydBqtYbjzz77LGQyGXbt2oVdu3aVuU5ERAR+/fVXNGrUCBMnTrT8CyMiIiKbZvVzopycnLB48WJMnToV8+fPR2RkJHx9fXHmzBlkZGQgICAAs2fPLnNOeno6hg0bBgA4cOAAfH19AQBdu3bFrFmzsGzZMrz11lv47rvv0KpVK1y5cgVXr16Fs7MzVq9ejaZNm9b46yQiIiLbYvU9UQAQHByMbdu2YciQIbh16xYOHjyIBg0a4OWXX8aGDRvg6uoquqypU6fihx9+QP/+/XHnzh3ExMRAo9Fg7Nix2L17N3r06GHBV0JERER1hdl6ojQaDa5duwaNRoOGDRuiVatWZp19r1KpsHLlSlF5fX19kZiYWGF6SEgIQkJCzFU1IiIiqoeqHUQdO3YMX3/9NeLj48ss5e7o6Ih+/frh5ZdfRocOHap7GSIiIiKrUq3hvFWrVmHy5Mk4cuQIioqKIAiC4Uej0WD//v0YO3Ys9u7da676EhEREVkFyT1RR44cwapVqwAAI0aMwNixY+Hv7w8XFxfcv38fFy9eREREBGJjY/H222/Dz88PKpXKbBUnIiIiqk2Se6IiIiIgk8kwc+ZM/Pe//0WPHj3QoEED2NnZwc3NDSEhIfjqq68wefJkaLVarFu3zpz1JiIiIqpVkoOoM2fOoEmTJpg6dWql+WbNmgU3NzccO3ZM6qWIiIiIrI7kIKqgoADe3t5V7i+jUCjQokUL5OXlSb0UERERkdWRHES1a9cOKSkp0Gg0lebTarVITU1Fy5YtpV6KiIiIyOpIDqImTZqEvLw8LFy4sNJ8q1atQm5uLp599lmplyIiIiKyOqKezjt69Gi5Y82aNUPPnj2xbds2XL16Fc8++ywCAgLg4uICtVqNq1evIjIyEgcOHMCjjz6KsLAws1eeiIiIqLaICqImT55c6dynU6dO4dSpUxWmR0dHY//+/bh48aLpNSQiIiKyQqLXiRIEoVoXqu75RERERNZEVBB16dIlS9eDiIiIyKZUa9sXIiIiovqKQRQRERGRBKKG8+bMmQOZTIa33noL7u7uhmOmkMlkWLJkiek1JCIiIrJCooKoXbt2QSaT4aWXXjIEUfpjVU0Y1+dhEEVERER1iaggauTIkZDJZGjQoEG5Y0RERET1kaggavHixaKOEREREdUXnFhOREREJAGDKCIiIiIJRA3nTZw4sdoXkslk+PHHH6tdDhEREZE1EBVExcfHV/okntin9IiIiIjqCpOezqvIzp074e7ujoceeshsFSOi2qPTFUOrK8aDfC3y1Fq4Oivg4qSAwt4OcjlnARARAdV4Oq+0nTt3olWrVli0aJFZKkVEtUdTWIT4C3cQGZuEq2m5huNtfdwwZqAfegV6wclB9N7lRER1Ft8JichAU1iE9VEX8NuRa+XSrqblYunGUxgW2hqTHw9kIEVE9R775YkIQMkQXvyFO0YDqNL2HrmGExfSodMV10zFiIisFIMoIgIAaHXF2B5zRVTeyNgr0DKIIqJ6jkEUEQEAHuRrkXLrnqi8yWm5eJCvtXCNiIisG4MoIgIA5KlNC4ryGEQRUT3HIIqIAACuzgrT8itNy09EVNcwiCIiAICLUoG2Pm6i8rbzcYOLE4MoIqrfRD2jfOvWrSrzFBYWVpnP29tbXK2IqMYp5HYYM9APSzeeqjLvmIHtobDndzAiqt9EBVGDBw+uNF0mk+H8+fOV5pPJZLh48aJptSOiGiOX26FXoBeGhraudJmDYaGt0TPQkyuXE1G9J+pdUBAEs/xUR0pKCmbPno2BAweic+fOeOSRR/D5559DrVZXq1wAWLJkCfz9/fHll19WuywiW+bkYI/nHw/EnPE9yg3ttfVxw5zxPbjQJhHR/yfqnTAiIsLS9ajU2bNnMWnSJKjVanTu3BlBQUE4ffo01q5di9jYWGzevBmurq6Syv7jjz/w/fffm7nGRLbLycEeoZ2bo2egZ8neeflauCq5dx4R0b+JCqJ69epl6XpUqKioCLNmzYJarcaCBQsQHh4OANBoNJg5cyZiYmKwfPlyfPDBByaXfffuXbz99tvV7iUjqowtbuYrl5fUzcnBHk3dlLVdHSIiq2T1ffJ79uxBamoqQkJCDAEUADg5OWHhwoUYNGgQtm7dihkzZsDNTdyTRXrvvfcesrOz0a1bN5w+fdrcVSfiZr5ERHWYSe/eWq0WBw8eREJCAh48eIDmzZujf//+6NChg6Xqh5iYGABAWFhYubTGjRsjODgYsbGxOHz4MIYPHy663E2bNiE2NhYzZsxAdnY2gygyO27mS0RUt4keSzh79iwee+wxvP7661i/fj1+/vlnrFixAqNGjcKbb75plgnexly+fBkA4O/vbzTdz88PAHDp0iXRZV65cgVLlixBt27dMG3atOpXkuhfuJkvEVHdJ+rrb1ZWFl588UXcu3cPgiDA0dERLi4uuHv3LgBg79690Gg0WL16tdkrmJGRAQDw9PQ0mu7h4VEmX1UKCgowa9YsKBQKLF26FHK53DwVrYQgCNDpdBWml06rLB9JV9NtrC0STNrMt2eAJwDbnZvHe9iy2L6Wxfa1PFtrY7FzpUU/nZebmwsvLy98/PHH6N+/PwAgLy8P69evx9dff42YmBicO3cOQUFB0mtthL6Hy8nJyWi6/rjYnrD//ve/uHz5MpYsWQJfX1/zVLIK+fn5SEhIEJX33Llzlq0MWbyNlUol3L1ambSZ770HBfj7znXk5+dbtG41gfewZbF9LYvta3l1qY1FDefFxcXBzs4OX3/9tSGAAgBXV1e8/vrrmDJlCgRBwKFDh8xeQX1PkUwmqzSfmKjx4MGD2LhxI4YNG4aRI0eao3pE5cjlctxXF5p0Tl5+YY30ihIRkfmI6om6efMmWrRoAZVKZTR9xIgRWLduHa5cETd8YQoXFxfk5ORU+A1do9EAKPn2X5m///4b7777Lpo3b46PP/7Y7PWsjFKprHBOF1DStamPzIOCgvhhagE13cY5eaYFUW6uTmjk5Weh2lge72HLYvtaFtvX8mytjRMTE0WNDIgKoh48eIDWrVtXmN6iRQsAwP3798XVzgQeHh7IyclBZmam4Tql6edC6edGVWTNmjW4e/cuOnbsiPnz55dJu3DhAgBg//79uH79Otq1a4fp06eb6RWU9KKJvWHkcrnV31y2ribaWL+Zb+llDSqi38y3rvzeeQ9bFtvXsti+lmcLbVzV6JeeqCCqqKgI9vYVZ3VwcABQMmnb3Pz9/XH58mUkJSWhW7du5dKTkpIM+SqjnzP1119/4a+//jKa5/Lly7h8+TJ69epl1iCK6h9u5ktEVPdZ/Tu3fg5WdHR0ubTs7GwcP34cCoUCffr0qbScxYsXIzEx0ejPxIkTAQCvvvoqEhMTsWHDBvO/EKpXSm/mWxlu5ktEZLus/p07LCwM3t7eiIuLw6ZNmwzHNRoN5s6dC7VajfDwcLi7uxvStFotkpOTkZycDK1WWxvVJuJmvkREdZzVv3s7OTlh8eLFmDp1KubPn4/IyEj4+vrizJkzyMjIQEBAAGbPnl3mnPT0dAwbNgwAcODAgRpbyoDo37iZLxFR3SU6iMrKysIvv/xSrTxSlxUIDg7Gtm3bsGrVKsTHxyMpKQm+vr4IDw/HlClT4OrqKqlcoprAzXyJiOom0UHU9evX8e6771aYLpPJKs0jk8mqtTaTSqXCypUrReX19fVFYmKi6LLnzp2LuXPnSq0aERER1UOigyixS6Bb6nwiIiIiayIqiDJlc18iIiKi+qDGZrU+ePCgpi5FREREZHEWD6IuXryI999/H/369bP0pYiIiIhqjEWWOCgoKMCePXuwZcsWnD9/3hKXICIbo9MVQ6srLlnqQa2FqzOXeiAi22bWICo5ORk//fQTfv31V8M+eoIgQC6XY+DAgea8FBHZEE1hEeIv3EFkbFKZ/QTb+rhhzEA/9Ar04qKjRGRzqv2uVVRUhP3792PLli04efIkgH+exGvWrBmefPJJPPXUU/D09KzupYjIBmkKi7A+6gJ+O3KtXNrVtFws3XgKw0Jbc/V2IrI5kt+xbt68iZ9//hk7duzA3bt3AfwTPLm4uGDhwoV4+OGHrX6nZiKyHJ2uGPEX7hgNoErbe+QaOrV1R2jn5hzaIyKbYVIQJQgCYmJi8NNPP+GPP/6AIAiGwKlVq1YYMWIEvvzySzRo0ABDhgyxSIWJyHZodcXYHnNFVN7I2CvcjJmIbIqoICojIwNbt27F9u3bkZ6ebgicmjRpgmHDhmHEiBHo3LkzAODLL7+0XG2JyKY8yNci5dY9UXmT03LxIF/LIT0ishmi3q0GDhyI4uJiCIIAFxcXDB48GMOHD0efPn04XEdEFcpTa03Ln6/l/oJEZDNEBVE6nQ4ymQwDBgzA9OnT8Z///MfS9SKiOsDVWWFafqVp+YmIapOoyQetW7eGIAg4dOgQnn76afTr1w/z58/HqVOnLF0/IrJhLkoF2vq4icrbzscNLk4MoojIdogKovbt24fNmzdj1KhRUCqVyMjIwJYtWzB+/HgMHjwYX3zxBVJSUixdVyKyMQq5HcYM9BOVd8zA9lDYc1I5EdkO0e9Y3bp1w8KFC/HHH39g4cKF6NatGwRBQFpaGtauXYthw4ZhzJgxlqwrEdkYudwOvQK9MDS0daX5hoW25pN5RGRzTH7HUiqVGD16NDZt2oT9+/dj6tSp8PDwgCAIuHDhAmQyGTIzMzF9+nT8/vvv0Ol0lqg3EdkIJwd7PP94IOaM71FuaK+tjxvmjO/BhTaJyCZV612rZcuWmDVrFt544w3ExcVh+/btiI2NhVarxcGDB3Hw4EE0adIEI0aMwJgxY+DnJ65bn4jqFicHe4R2bo6egZ4le+fla+Gq5N55RGTbzPLVz87ODv369UO/fv2Qk5ODX3/9FTt37sSlS5eQlZWFH374AT/++CMuXrxojssRkQ2Sy0uCJScHey5jQER1gqivf4MHD8bMmTNFFdioUSNMmjQJv/zyC3bs2IFx48ahQYMGhgU6iYiIiOoCUT1RaWlp8PLyMrnwgIAABAQE4O2338bvv/9u8vlERERE1qpGJiI4ODhg2LBhNXEpIiIiohrB2ZxEREREEjCIIiIiIpKAQRQRERGRBKKXODh9+jQ6duwo+UIymYxLHBAREVGdYdI6UVymgIiIiKiE6CCqVatWmDZtmiXrQkRERGQzRAdRTZs2xahRoyxZFyIiIiKbwYnlRERERBIwiCIiIiKSgEEUERERkQQMooiIiIgkEDWxPCIiAg0aNLB0XYiIiIhshqggqlevXpauR5VSUlKwevVqnDp1CllZWfDy8sLQoUMxbdo0ODs7m1TW0aNH8d133+HPP/9EQUEBmjdvjocffhgvvvgiGjVqZJkXQERERHWKSYtt6hUWFuLkyZM4fvw4bt++jZycHMhkMjRs2BDt2rVDt27d0LNnT8hkMrNU8uzZs5g0aRLUajU6d+6MoKAgnD59GmvXrkVsbCw2b94MV1dXUWVFRkbivffeAwB06dIF7u7uuHjxIr799lvs2rULW7Zsga+vr1nqTURERHWXSUGUVqtFREQEvvnmG+Tm5laa18PDA1OnTsUzzzwDOzvpU6+Kioowa9YsqNVqLFiwAOHh4QAAjUaDmTNnIiYmBsuXL8cHH3xQZVmpqan48MMPoVAo8PXXXyM0NNTwuj788ENERkZi/vz5+PrrryXXl4iIiOoH0dFNXl4epkyZgs8++ww5OTkQBAEuLi5QqVTo2rUrOnXqhFatWkEul0MQBKSnp+PTTz/Fiy++CLVaLbmCe/bsQWpqKkJCQgwBFAA4OTlh4cKFcHZ2xtatW6sM6gAgKioKWq0WY8aMMQRQAKBQKPDmm28CAOLi4lBYWCi5vkRERFQ/iO6JmjFjBuLj4yGXy/HUU0/hqaeegr+/f7l8hYWFOHv2LLZt24aoqCgcOXIEc+bMwapVqyRVMCYmBgAQFhZWLq1x48YIDg5GbGwsDh8+jOHDh1da1ksvvYThw4cbnUOl0+kAAHZ2dtXqOSOqy3S6Ymh1xXiQr0WeWgtXZwVcnBSQy2WQy+WGvyMiovpAVBAVGxuLP/74A66urli7di169OhRYV4HBwf06NEDPXr0wJgxY/DSSy/hwIEDOHbsGHr37m1yBS9fvgwARgM2APDz80NsbCwuXbpUZRBlZ2eHli1bljuenZ2Njz76CADwxBNPwN5e0lQxMqOKPqwV9naQyxnk1gZNYRHiL9xBZGwSrqb90/Pb1scNYwb6oavKH9dTkmuxhkRENUtUtLBr1y7IZDK8++67lQZQ/9arVy+8+eab+OSTTxAVFSUpiMrIyAAAeHp6Gk338PAok88UX3/9NeLi4pCQkIDCwkKMHDlS1NwqUwmCUOk39NJp/CYPaIsExF+s+MO6V4AXFPamPbTANq4erU7A+qgL+O3ItXJpV9NysXTjKQwNbY0Jj6rYvhbA+9ey2L6WZ2ttLAiCqHyigqiLFy/CwcEBI0aMMLkio0aNwsKFC3H27FmTzwVgmE/l5ORkNF1/XMq8q7179+Kvv/4CAMjlcty7dw83b95Eu3btJNW1Ivn5+UhISBCV99y5c1XmUSqVhqGT/Pz8atbOesjlcvi174Aff7uEfUevl0v/94d18pVESX+MYtqY/uHh4YXkDJ3RAKq0345cQ2CbJmjbTIaMjPSaqVw9xPvXsti+lleX2ljUuEhmZiZatGgBhUJh8gWcnZ3h6+uL27dvm3wuUPLBCqDK5RLERo2lrV27FmfPnkVUVBQeeeQRxMTE4Nlnn5VcV0uSy+Xw8PCCf8dAuHu1gr2LB9y9WsG/QyA8PLwM7WTLmjZthpOXMowGUKX9duQaTif+jaZN3WuoZvVb46ZNseOguGG6nQeT0Zi/FyKqJ0T1RBUUFIheh8kYNzc33Lx5U9K5Li4uyMnJqbDHRaPRACjpnTGVl5cXAEClUmHFihW4f/8+4uLisH79esydO1dSfY1RKpUVzukCSro29ZF5UFCQ0YBIP8T12bZjFQxxdTJ5iMvaaIsELN16WFTenQeTsfiVvvD2bi4qv5g2JuNy8gqRcuueqLzJabkoKBTQpUsXy1aqnuH9a1lsX8uztTZOTEwUNdIjKojS6XTVesH29vYoLi6WdK6HhwdycnIMvWH/pp8LpZ8bVR0jRoxAXFwcLl68WO2ySpPJZKLbTy6Xl8urKSzC+t2Vz0cZFtoakx8PhJOD7U6Kz8nLN+nD+oFGi6ZupgfPxtqYKpan1pqWX6NF00am/15IHN6/lsX2tTxbaGOxi4Vb/WNO+h6cpKQko+n645X19OhFRETgzTffxPnz542mOzg4AChZ4NNa6HTFiL9wp8r5KHuPXMOJC+nQ6aQFq9bA5A/rfNPykzSuzqYN47sqTR/2JyKyRVYfRPXv3x8AEB0dXS4tOzsbx48fh0KhQJ8+faosKyEhAbt378aOHTuMph86dAhASVejtdDqirE95oqovJGxV6C14SCKH9bWyUWpQFsfN1F52/m4wcWJvxciqh9Ej/3cv38fJ06ckHSR+/fvSzoPKFlk09vbG3Fxcdi0aRPGjRsHoGQu1Ny5c6FWq/HMM8/A3f2fyaxarRY3btwAALRs2dIwIf7ZZ5/Fnj178NNPP6Fv374YNGgQgJJJ6Rs3bsTOnTvh6OiIiRMnSq6vuT3I15o2xJWvtdkhPf2Hdek5XxXhh3XNUcjtMGagH5ZuPFVl3lED2kEut+25eUREYon+tL1y5UqtBBdOTk5YvHgxpk6divnz5yMyMhK+vr44c+YMMjIyEBAQgNmzZ5c5Jz09HcOGDQMAHDhwwLChcI8ePTBz5kx8/vnnmD59Ojp16gQvLy9cvnwZN27cgFKpxLJly4wuyFlbpAxxSZknZA1M+bAeM7A9FPZW35FaJ8jldugV6IWhoa0rHVYeGtoa3fzdYccYiojqCdGfQoIgVOunOoKDg7Ft2zYMGTIEt27dwsGDB9GgQQO8/PLL2LBhg0lPDr700kv4/vvv0a9fP6SmpuLQoUMoKipCeHg4fvnlFwwePLhadTW3+jTEVfrDujLDQlujZ6AnVy6vQU4O9nj+8UDMGd+j3NBeWx83vDW+OyY8quKK5URUr4jqiTpw4ICl61EllUqFlStXisrr6+uLxMTECtNDQ0PLbEBszerbEJf+wzqorTu2x14pt5xD+MD26BnoabNDlrbMycEeoZ2bo2egZ8l2PPlauCr/2TvvwvlzNrESMRGRuYj6JPLx8bF0PagC9XGIq7IPa+6dV7vk8pL2d3KwLzNsrNPpGEARUb3Dr/NWTux8lLo2xFXRh7Wt4obKRER1D4MoG8AhLtumKSxC/IVKNlQO9OLvjojIBvGd20ZwiMs2aQqLsD6q7q82T0RUH/Fd24bUtSGuus6U1eY7tXVHaOfmDIaJiGwI37HJauh0xdAUFiErNx/Xb99DVm4+NAVFNruVTX1abZ6IqD5iTxRZhbo4b6g+rTZPRFQfsSeKap1+3tDSjafKrYelnzf0fdQFaAqtZ2NoMbihMhFR3cYgimqVKfOGTlxIt6mhvfq02jwRUX1U7bGD4uJipKSk4N69e9DpdJVu8dKzZ8/qXo7qGFPnDdnSWlj1bbV5IqL6plpB1M8//4wVK1YgJyenyrwymQwXL16szuWoDqrL84bq42rzRET1ieR37QMHDuDDDz9Edna2qA2Ii4ttZxiGak5dnjfEDZWJiOo2yV/pN2zYAAAICQnBW2+9hXbt2sHR0dFsFaP6oa7PG+Jq80REdZfkd+4LFy7A2dkZX375JVxdXc1ZJ6pH6sO8Ia42T0RUN0l+99ZqtWjTpg0DKKoW/bwhMWx53lDpleZbeTVEUzclnBztGUAREdkwye/grVq1wu3bt81ZF6qHOG+IiIhsleRPpCeeeAJ3797F3r17zVkfqof084bmjO+Btj5uZdLa+rhhzvge3KCXiIisjuRPpYkTJ+Lw4cN4//33kZmZiYEDB8LT0xMKRcVzVuzs2ItAxnHeEBER2RrJQdT48eNRWFiIBw8eYPHixVi8eHGl+blOFFVFLrcrM3eIiIjImkkOohISEgz/X9kq5URERER1keQgKiIiwpz1ICIiIrIpkoOoXr16mbMeRERERDaFs3WJiIiIJKj2M+PJyclITEyERqMptz+eTqdDfn4+0tPTcejQIS6HQCSSTlcMra645ElFtRauznxSkYjI2kgOooqLi/H2229j9+7dVeYVBAEymUzqpYjqFU1hEeIv3EFkbFK5vfbGDPRDr0AvrplFRGQFJL8TR0ZGIioqCgCgUCjg5uaGv//+G25ubnB0dER2dja0Wi1kMhk6duyICRMmmK3SRHWVprAI66Mu4Lcj18qlXU3LxdKNpzAstDUXHyUisgKSxwV2794NmUyGiRMnIiEhAdHR0bC3t0f//v3xf//3fzh16hTmz58PBwcHZGRkYMCAAWasNlHdo9MVI/7CHaMBVGl7j1zDiQvp0OmKK81HRESWJTmISkxMhFKpxKxZsyCXy+Hi4gKVSoXjx48DABwcHDB27Fi89dZbyMrKwoYNG8xWaaK6SKsrxvaYK6LyRsZegZZBFBFRrZIcROXl5cHX1xdOTk6GY35+fsjIyMDdu3cNx8aOHQtnZ2ccPHiwWhUl6XS6YmgKi5CVm4/rt+8hKzcfmoIi9mRYmQf5WqTcuicqb3JaLh7kay1cIyIiqozkSRVKpbLcXngtW7YEAFy9ehVNmjQBUNIj1apVK6SmplajmiQVJynbjjy1aUFRXr6W2+MQEdUiyT1R3t7eSE1NRUFBgeGYr68vBEFAYmJimbyFhYUoLCyUXkuSRD9JeenGU2UCKOCfScrfR12AprColmpIpbk6V7x5t9H8StPyExGReUkOonr16gW1Wo3//ve/hvWhOnbsCKBk0rl+P73k5GRcu3YNXl5eZqguicVJyrbHRalAWx83UXnb+bjBxYlBFBFRbZIcRE2YMAEKhQKbN2/GgAEDUFhYCH9/f3Ts2BEJCQmYMmUKlixZgueeew7FxcXo2rWrOetNVeAkZdujkNthzEA/UXnHDGwPhT0X3SQiqk2S34VbtmyJ5cuXo0GDBnjw4AEcHBwAAG+88QZkMhmOHj2KH374AZmZmXBycsIrr7xSrYqmpKRg9uzZGDhwIDp37oxHHnkEn3/+OdRqtcllHTx4EC+88AJ69+6NwMBA9OnTB6+99hr+/PPPatXRmnCSsu2Ry+3QK9ALQ0NbV5pvWGhr9Az05MrlRES1rFrvwg8//DB+//13LF261HCsf//+WL9+Pfr06YNWrVph0KBB2Lhxo2HSuRRnz57F6NGjERUVBXd3dwwYMABqtRpr167F008/jby8PNFlLV++HNOmTUNcXBx8fHwwYMAANGzYEPv378czzzyDnTt3Sq6nNZEySZlqn5ODPZ5/PBBzxvcoN7TX1scNc8b34EKbRERWotrvxA0bNsSgQYPKHOvduzd69+5d3aIBAEVFRZg1axbUajUWLFiA8PBwAIBGo8HMmTMRExOD5cuX44MPPqiyrJMnT2LdunVQKpVYt24dgoODDWk//fQTPvzwQ3zwwQcIDg6Gt7e3WepfWzhJ2XY5OdgjtHNz9Az0LNk7L18LVyX3ziMisjZW/268Z88epKamIiQkxBBAAYCTkxMWLlwIZ2dnbN26Fbm5uZWUUmL79u0AgBdeeKFMAAUATz/9NPr374/CwkJER0eb90XUAk5Stm1yuR2cHOzR1E2JVl4N0dRNCSdHewZQRERWRFRPlD74qK7SQZBYMTExAICwsLByaY0bN0ZwcDBiY2Nx+PBhDB8+vNKynJycoFKpKuwla9u2LQ4dOoT09HST62lt9JOUl248VWVeTlImIiIynaggat68eZDJZNW+mJQg6vLlywAAf39/o+l+fn6IjY3FpUuXqgyiPvroo0rT9RPLzb0cgyAI0Ol0FaaXTqssn6l6BZRMUq5smYNhoa3RM8ATQOV1tHWWamMqwfa1LLavZbF9Lc/W2li/TFNVRAVRFc0P0mg0hi1eGjRoAD8/P7i5uUGj0SApKQl///03ZDIZOnTogEaNGomr+b9kZGQAADw9PY2me3h4lMknVUxMDE6fPg2FQmG016s68vPzkZCQICrvuXPnzHZdR0dHTHhUhcA2TbDjYHK5FctHD2iHbv7uSLpyqcyiqXWdOduYymP7Whbb17LYvpZXl9pYVBClH1IrLS8vD8888wwKCwsxb948DB8+HPb29uXOmzdvHvLy8vDdd99JqqB+CYPSe/SVpj8uZakDvcTERLz77rsASuZL+fj4SC7LmhQUFCD5SiLaNnPHgpd6Q63RIS+/EK5KBzg7ypF9928kX0m0iW8FFVEqlZDL5dDpdMjPz6/t6hARUT0i+em81atXIykpCV9//TUeeugho3kGDRqElStXYvz48fj888/xySefmHwduVyO4uLiKocTxXa9/dvZs2cxdepU5OTkYODAgXj99dcllVMZpVJZ4XAkUNK1qY/Mg4KCIJfLzV4HAFA6KtDU7Z9gVOndHN7ezS1yLUsrFgCdTsADjRb31Vo0cC55ek0ul8HOyK1SU21cX7F9LYvta1lsX8uztTZOTEwU9cVcchC1b98++Pr6VhhA6fXo0QOtW7dGbGyspOu4uLggJyenwhej0WgAlAQqptq3bx/eeecd5Ofn4+GHH8bnn39eblNlc5DJZKJvGLlcbvU3V22r7qbKbGPLYvtaFtvXsti+lmcLbSx2HrjkiOHu3bto2LChqLxKpRIPHjyQdB39nKfMzEyj6fq5UPp8Yq1evRpvvPEG8vPz8eyzz2LlypWGVdfJenFTZSIishaSgygvLy9cuXIF2dnZleZLTU3F5cuX4evrK+k6+mGwpKQko+n645UNl5VWXFyMd955BytXroSdnR3mzp2LDz/80OqjYuKmykREZF0kB1GDBg1CYWEh3nzzTdy/f99onoyMDMyYMQPFxcUYNmyYpOv0798fAIwugJmdnY3jx49DoVCgT58+osqbN28edu7cCWdnZ3z11VeYOHGipHpRzeOmykREZE0kz4l6/vnnERUVhaNHj+Lhhx/G4MGDoVKp4OzsjLy8PFy8eBExMTFQq9Xw8/PDhAkTJF0nLCwM3t7eiIuLw6ZNmzBu3DgAJXOh5s6dC7VajWeeeQbu7u6Gc7RaLW7cuAGgZKNkhaJkNe5ffvkFkZGRsLe3x5o1axASEiL15VMtkLKpMveYIyIiS5H8CdOsWTN8++23ePPNN5GcnFxu417903LdunXDsmXL4OrqKuk6Tk5OWLx4MaZOnYr58+cjMjISvr6+OHPmDDIyMhAQEIDZs2eXOSc9Pd3Q83XgwAH4+vpCp9NhxYoVAICmTZsiMjISkZGRRq/Zt29fjBw5UlJ9yXKkbKrc1M30Bw6IiIjEqNbX9A4dOmDXrl2Ijo7GoUOHkJKSgtzcXDRu3Bht2rTBI488ggEDBlR7tfPg4GBs27YNq1atQnx8PJKSkuDr64vw8HBMmTJFVICWmJiI27dvAygJsqKioirM27BhQwZRVoibKtddOl0xtLrikg2X1Vq4OnPD5X+Ty+Vo2rQZtEUCcvLy2U5EVqDaYx1yuRzDhg2TPOdJLJVKhZUrV4rK6+vri8TExDLHAgICyh0j26LfVPnfT+UZw02VbUd1l6yoL9q198epS5n4bFsc24nISvCrC9kM/abKYnBTZdvAJSvE0eoERPyWiGWbz7CdiKyIqK8tX3zxBQBg0qRJhj3w9MdMMWPGDJPPIdKTy+3QK1DkpsqBnhzesHKmLFnRqa07Qjs3r5e/U3077Tt6vdJ89b2diGqDqCDqq6++gkwmw4gRIwxBlP6YGIIgQCaTMYiianNysMfzjwciqK07tsdeKTesET6wPXoGenJYwwaYumRFfQ2M2U5E1kvUJ03Pnj0BlN1aRX+MqKY5OdgjtHNz9Az0LJmInK+Fq5ITbG0Nl6wQh+1EZL1E/aVt2LBB1DGimiKXlwRLTg72XMbARnHJCnHYTkTWi1/ZiahWcMkKcdhORNbLLEFUbm4u9u3bhzVr1mDJkiUAgIKCAsTHx5ujeCKqg/RLVohRn5esYDsRWa9qBVGCIGDlypUYMGAAZs6ciS+//BI//PADAODmzZuYNGkSnnnmGdy9e9ccdSWiOoRLVojDdqqYTlcMTWERsnLzcf32PWTl5kNTUMTNx6nGVGv24Zw5c7B7924IgoDGjRujoKAA+fn5AICcnBwIgoCEhARMmDAB27dvLzMxnYjqNy5ZIY5cbodeAV54NKRVpcsc1Ld24iKtZA0k/7Xt378fUVFRaNKkCb755hscPXoUHTp0MKR3794dmzZtQtOmTXH16lVERESYpcJEVHfol6yYM75HuSGrtj5umDO+ByY/HljvPwwV9jJMHOqP2eO6sZ3ARVrJekj+i/v5558hk8mwbNky9O7d22ie7t27Y8WKFRg/fjyio6Mxbdo0yRUlorqJS1aIk3wlEW2buWPJK33xQFN/24mLtJI1kRxEnT9/Hs2bN68wgNLr0aMHfHx8cO3aNamXIqI6jktWVE2n0yEjIx3e3s3R1E1Zb9uJi4+SNZF8Z6nVasPq5VVp0qQJiorYrUpERNUjZfFRIkuR3BPl7u6O69evG7Z0qYhWq8W1a9fg7u4u9VJERJXS6Yqh1RWXDAeqtXB1rn/DXPUFFx8layI5iOrVqxd27dqFzZs3Y9y4cRXm+/HHH3H//n0MHDhQ6qWIiCrEp7TqFy4+StZE8le0yZMnw87ODkuWLEFERASys7PLpGdlZWHFihVYvnw57OzsMH78+GpXloioND6lVf9w8VGyJpKDqA4dOuC9996DVqvFokWLEBoaijNnzgAAQkJC0LdvX6xbtw7FxcV47bXX0LlzZ7NVmojIlKe0TlxI5wKMdQQXHyVrUq27a9y4cVi3bh38/f0hCILhJzs7G4IgoGXLlli+fDmmT59urvoSEQEw/SktLYOoOqH0Iq2VqW+Lj1LtqPZEgX79+qFfv35IS0vDlStXcP/+fSiVSrRp0wbt2rUzRx2JiMqR8pSWJeZGcVJ7zdMv0hrU1h3bY6+UmwsXPrA9egZ6ci4cWZzZ7jAfHx/4+PiYqzgiokpZw1NanNRee7hIK1kD/nUTkU2q7ae09JPajc3J0k9qHxbaul5tx1LTuEgr1TZRf9kDBgyo9oVkMhliY2OrXQ4REfDPU1r/firPGHM/pcWtR6wTh1appokKou7cuQOZTAZBECRfqLIFOYmITKV/SmvpxlNV5jX3U1rcesT6cGiVaoNJd5RMJkPHjh0xdOhQrkBORLWq9FNalfUIWeIpLWuZ1E4lOLRKtUXU3fTWW2/ht99+w/nz53Hx4kUkJiaiV69eGDZsGIYMGYKGDRtaup5EROXU1lNa1jCpnUpwaJVqk6h3lilTpmDKlClITU3F3r178dtvv+Ho0aM4duwY5s+fjz59+uCxxx7D4MGD4ezsbOk6k5lw/gDVBbXxlFZtT2qnf3BolWqTSV/PWrRogWnTpmHatGlISUnBnj178Ntvv+HgwYM4dOgQHB0dMWDAADz22GPo378/HBwcLFVvqibOH6C6pKaf0qrNSe1UFodWqTZJDsfbtGmDV199FXv27MGuXbswdepUeHp6Yt++fXj99dcREhKCt99+G4cOHYJOpzNnnamauN8YUfVw6xHrIWVolchczPKXrVKpMHPmTERHR2PHjh14/vnn4ebmhl9//RUvvfQSQkND8cEHH5jjUlRN3G+MqPq49Yj14NAq1Saz/2UHBATgrbfeQkxMDN5//304OzsjNzcX27ZtM/elSALuN0ZkHvpJ7XPG90BbH7cyaW193DBnfA8+DVYD9EOrYnBolczN7H/df/75J/bt24fo6Gjcvn3bsLYUt4SxDpw/QGQ+3Hqk9tXmemFEZvl0TEhIMAROd+7cMQRO3t7eGDJkCIYOHYrOnTub41JUTXw0m8i8uPVI7arN9cKIJAdRZ86cwb59+7B///4ygVPz5s0NgdN//vMfs1WUzKOi+QN2djK09GwAZyd7qDVFuJF+H8XFAucPEJHVq631wohMuqNOnz5tCJzS09MNgZOnp6chcOratatFKpqSkoLVq1fj1KlTyMrKgpeXF4YOHYpp06ZVa22q7OxsPP744wgMDMS6devMWGPr9O9Hs12VCgwLbYNHQ1qhQKvDvQeFaOjiAAd7OeIv3qnW/AGuQ0VENYVDq1QbRAVRCxYswP79+5GRkWEInJo1a4YhQ4bg0UcfRY8ePSxaybNnz2LSpElQq9Xo3LkzgoKCcPr0aaxduxaxsbHYvHkzXF1dTS5XrVbjtddeQ2ZmpgVqbZ1Kzx9o7u6Cj18MwYWrWfj0+/hy395G9W8n+Tpch4qIahqHVqmmifoU27BhA2QyGeRyOXr37o2hQ4eiR48esLMriexTU1NFXaxFixYmV7CoqAizZs2CWq3GggULEB4eDgDQaDSYOXMmYmJisHz5cpOXUEhNTcXMmTNx7tw5k+tky/TzB0b2b4fH+rRBZMwV7Dt2vVy+q2m5WLb5tKT9priPFRER1QcmfYLpdDr88ccf+OOPP0y+kEwmw8WLF00+b8+ePUhNTUVISIghgAIAJycnLFy4EIMGDcLWrVsxY8YMuLlV/ZirRqNBREQE1q1bh7y8PLRo0UJ0EFhXODnYY/yjHXD03G2jAVRppu43xX2sbAuHXImIpBMdROmH8aSSen5MTAwAICwsrFxa48aNERwcjNjYWBw+fBjDhw+vsry9e/di2bJlaNSoERYsWAA7Ozu8++67kupmywQAOw4micpryn5T3MfKdnDIlYioekS9Qx44cMDS9ajQ5cuXAQD+/v5G0/38/BAbG4tLly6JCqIaNWqEl19+GZMnT0bDhg2xY8cOs9bXGEEQKt36pnRaTW2RI2W9KIVcVmvlVldttLE10+oEcUOuwwOhsK/698P2tSy2r2WxfS3P1tpYbMePqCCqNhfKzMjIAFDyBKAxHh4eZfJVZdCgQRg0aJB5KidSfn4+EhISROWtiTlarq6usHfxMOmc3DwN/r5zA3l5eTVerrnVt3lw/+bh4YXkDJ2oIdeANk3QtpkMGRnposuv7+1raWxfy2L7Wl5damOrH0dRq9UASuZAGaM/rs9HVdPpdGjg7GDSOa5Khyq/PViqXDKvxk2bYsfBZFF5dx5MRuOm7hauERGRbbL6CQ9yuRzFxcWQySofUqjunC1LUiqVFQ5HAiXBhz4yDwoKglwut3idtEVCmfWiKtPOxw0NXRzRtJLXYOlyq6s22tha5eQVmjTkWqAV0KVLl0rzsX0ti+1rWWxfy7O1Nk5MTER+fn6V+aw+iHJxcUFOTk6FL0aj0QAoCVSslX55CDHkcnkN3VzFJu83JW4CuKXKNZ+aa2PrZOmtf+p7+1oa29ey2L6WZwttXFXHjZ7VD+fp5zxVtCCmfi6UPh+JU3q/qcqYut+Upcol86lo658K83PrHyIio6z+E0w/DJaUZPxxfP3xyobLyDj9flNzxvdAW5+ya2y19XHDnPE9JC2IaalyyTz0W/+I0c7HrVpb/xAR1WVW/ynWv39/REVFITo6GmPHji2Tlp2djePHj0OhUKBPnz61VEPbZqn9priPlfUqvfVPVfRDrkREVJ7VvzuGhYXB29sbcXFx2LRpk+G4RqPB3LlzoVarER4eDnf3f54g0mq1SE5ORnJyMrRa0+Z/1Eel95pq5dUQTd2UcHK0r3agY6lyqXo45EpEZB5W3xPl5OSExYsXY+rUqZg/fz4iIyPh6+uLM2fOICMjAwEBAZg9e3aZc9LT0zFs2DAAJQuF+vr61kbViayWfsg1qK07tsdeKbdiefjA9ugZ6MkhVyKiStjEO2RwcDC2bduGVatWIT4+HklJSfD19UV4eDimTJkCV1fX2q4ikc3hkCsRUfXYRBAFACqVCitXrhSV19fXF4mJiaLyjh49GqNHj65O1YhsllxuV2bYlYiIxONXTSIiIiIJGEQRERERScAgioiIiEgCm5kTRUSWodMVQ6srLplcrtbC1ZmTy4mIxGAQRVSPaQqLEH/hDiJjk8otczBmoB96BXpxmQMiogrw3ZGontIUFmF91AX8duRaubSrablYuvEUhoW25hY9REQVYF89UT2k0xUj/sIdowFUaXuPXMOJC+nQ6YprpmJERDaEQRRRPaTVFWN7zBVReSNjr0DLIIqIqBz20RPVQw/ytUi5dU9U3uS0XDzI13JIjyTjwwtUV/FdkageylObtjF3Xr6WK5qTJHx4geoy3rlE9ZCrs8K0/ErT8hMBfHiB6j72oxLVQy5KBdr6uInK287HDS5ODKLINHx4geoDBlFE9ZBCbocxA/1E5R0zsD0U9nyrqI90umJoCouQlZuP67fvISs3H5qCIlEBDx9eoPqA/adE9ZBcbodegV4YGtq60p6CYaGt0TPQk5N/66HqzmXiwwtUH/COJaqnnBzs8fzjgQhq647tsVfKfVCGD2yPnoGeJn2wyeVyFAuAtrCIT2LZMHPMZeLDC1QfMIgiMlFdelzbycEeoZ2bo2egZ8nrydfCVSnt9Tg6OqJVm3Y4eu42n8SyYabMZerU1h2hnZsbvU/48ALVB3w3IzJBXXxcWy4vCZacHOyr1RPQuo0ffvztEvYdvV4uraLei7oUkNYVps5lqmi4V//wQum/k4rw4QWyVbb1bk9Ui/i4dsWKBeDkpQyjAVRppXsvtP+/x6MuBaR1gbnmMukfXli68VSV5fDhhbrH2Bck/w6ByL6bVdtVMyu+QxGJYK4hjrpKpxOw42CyqLyRsVfQM8AT3+++gL0MSK2OueYy8eGF+quyHvtR/duimYcAubwWK2hGvGuJRODj2pV7oDGt9+Leg0JcTLlbaT6uH1Q7zDmXSf/wwpzxPcqtS9bWxw1zxvdgoFzH6Hvsl248VW4o92paLpZtPoPvd1+AprColmpoXrxzSbL6NJ+Fj2tX7r4JvRd2djI8yNcioE0TAMCN9PsoLhaM5q1szg1ZhrnnMpnz4QWybvWxx77+vMuTWdXFCdaV4ePalWsgovfCVanAsNA2eDSkFQBgYPcWePyhtnCwlyP62HXsPZKCvPyy7VwfA9LaZom5TOZ6eIGsm7keSrAlfGcik9XHCdZ8XLtyLk6V9140d3fBxy+G4MLVLHz6fXy5wPvxvm2x/I3++PCbo7j994My59a3gLS2cS4TSVUfe+x595NJ6ut+WNxrrnJyuQyj+rc1muaqVGD+1BDsiL2CL34+Y3SexBc/n8GOg0n4+MWQcgFofQtIrQHnMpEUUnrsbR3/Asgk9bG7FuDj2lWxkwHdOzTDoyGtyi1zMCy0Dc4nZ2HfscqXP9h39Bo6tGqMoaGtse1AyT1WHwNSa8G5TGSq+thjz78CMomU7tq6oPQQR2Xq8xDH9ZRkTBzqj7fGdzf0XtjZyfBoSCvsOixu+YOow1cxNKQ17GQl/66PAak1KT2PqZVXQzR1U8LJ0b5e3t9UtfrYY8+eKDJJfZ5gbYm95uqSgoICJF9JREhQEHoFeuFBvhYFWh1kMplJgXeBVoeWXg0R0KZJvQ1IiWxRfeyxr5/v9iRZfeyuLY1DHJXT6XSwkwEKB3tDMHn9trgASi8vX4sXRnSCf+vG9TYgJbJFNfVQgjUtr8N3KDIJ98Pi49qmMjXwdndTonEDx3ofkBLZoqp67EcPaIfgwOaSvyBZ2/I6DKLIJPWxu5aqx9TA21WpYABFZMMq6rF3tJch++7fUNjLJJVrjcvr8J2KTMIJ1mQqfeAtBgNvorrh3w8lNHJ1QOKlC8jISJdUnrUur8N3KzIZ15AhUzDwJqLqstb9S/kpR5JwgjWZgk82ElF1WOtq6DbzjpWSkoLVq1fj1KlTyMrKgpeXF4YOHYpp06bB2dnZpLLS09OxZs0aHDlyBHfu3IG7uzsGDRqEV155BU2aNLHQK6h7OMGaTMHAm4ikstbldWziXevs2bMYPXo0oqKi4O7ujgEDBkCtVmPt2rV4+umnkZeXJ7qs1NRUjBkzBj/99BOcnJwwcOBAyOVybNy4EaNGjcKdO3cs+EqI6jcu3khEUljr8jpW/85VVFSEWbNmQa1WY8GCBdi2bRtWrlyJ33//HYMGDUJiYiKWL18uurx33nkHmZmZeOWVVxAVFYWVK1ciOjoaTz/9NO7cuYMPP/zQgq+GiIiITGWtq6FbfRC1Z88epKamIiQkBOHh4YbjTk5OWLhwIZydnbF161bk5lb9+PTJkydx8uRJtG7dGq+++qrhuFwux7x58+Dt7Y2DBw8iKSnJIq+FiIiITGetT/lafRAVExMDAAgLCyuX1rhxYwQHB0Or1eLw4cOiyxo8eDDs7Mq+dIVCgUGDBgEADhw4UN1qExERkZlY61O+Vj+x/PLlywAAf39/o+l+fn6IjY3FpUuXMHz4cFFlqVSqCssCgEuXLkmtrlGCIECn01WYXjqtsnwkHdvYsti+lsX2tSy2r+WZo40VchmeHx6ITm2bVrxieYAXFHJZtX+PgiCIymf1QVRGRgYAwNPT02i6h4dHmXw1VZYp8vPzkZCQICrvuXPnzHptKo9tbFlsX8ti+1oW29fyqtPGcrkcbdzdseCl3lBrdMjLL4Sr0gHOjnJk3/0bl/46X6OBsNUHUWq1GkDJHChj9Mf1+cSUpVQaf+zRlLKIiIioZul0OmRkpCMjIx1KpRJyuRx/5+qQn59fK/Wx+iBKLpejuLgYMlnle+2I6XqTy+UAYJayTKFUKiscjgRKbgp9ZB4UFGSoJ5kP29iy2L6Wxfa1LLav5dlaGycmJooKzKw+iHJxcUFOTk6FL0aj0QCouHfp32UBMEtZppDJZKJvGLlcbvU3l61jG1sW29ey2L6Wxfa1PFto46o6W/Ss/uk8/TylzMxMo+n6+Uv6fDVVFhEREdVvVh9E6YfBKlq7SX+8suEyS5RFRERE9ZvVB1H9+/cHAERHR5dLy87OxvHjx6FQKNCnTx/RZf3vf/9DcXHZHZ61Wq1hfaiBAwdWt9pERERUx1l9EBUWFgZvb2/ExcVh06ZNhuMajQZz586FWq1GeHg43N3dDWlarRbJyclITk6GVvvPpoXdunVDUFAQkpOTsXz5csMEcp1OhwULFuD27dvo168fOnbsWHMvkIiIiGyS1U8sd3JywuLFizF16lTMnz8fkZGR8PX1xZkzZ5CRkYGAgADMnj27zDnp6ekYNmwYgJLVx319fQ1pixYtwvjx4/HNN9/gwIEDaN++Pf766y/cuHEDPj4++PTTT2v09REREZFtsvqeKAAIDg7Gtm3bMGTIENy6dQsHDx5EgwYN8PLLL2PDhg1wdXUVXVb79u2xY8cOjB49Gvfv30dsbCwAYMKECdi6dWuFC3ESERERlSYTzL0oEhkkJCRAp9NBJpNVumyCIAiGZReUSqXoRytJPLaxZbF9LYvta1lsX8uztTbOz8+HIAiQy+Xo0qVLhfmsfjjPluknrwuCIHoV9NpadbU+YRtbFtvXsti+lsX2tTxbauN/P4T2bwyiLEihUECr1cLOzg6Ojo61XR0iIiISoaCgAMXFxVAoFJXm43AeERERkQQ2MbGciIiIyNowiCIiIiKSgEEUERERkQQMooiIiIgkYBBFREREJAGDKCIiIiIJGEQRERERScAgioiIiEgCBlFEREREEjCIIiIiIpKAQRQRERGRBAyiiIiIiCRgEEVEREQkAYMoIiIiIgkYRBERERFJYF/bFbAVKSkpWL16NU6dOoWsrCx4eXlh6NChmDZtGpydnU0qKz09HWvWrMGRI0dw584duLu7Y9CgQXjllVfQpEkTo+ecO3cOa9aswfnz53Hv3j20aNECI0eOxKRJk6BQKMrlX758OdatW1dhHQYMGFBpek2r7fYt7eTJk5gwYQJmz56NKVOmVJhv7969iIiIwNWrV6HT6dChQwdMmjQJjzzyiEn1rQm21r68f8W3b1FREX766Sfs2rULSUlJKCwshKenJ/r164dp06bBy8vL6HVs6f4FbK+NeQ+Lb9/i4mJs27YN27Ztw5UrV2BnZ4d27dph5MiRePrpp2FvbzxUsYZ7WCYIglBjV7NRZ8+exaRJk6BWq9G5c2c0b94cp0+fRmZmJvz9/bF582a4urqKKis1NRXPPPMMMjMzoVKp0KZNG1y8eBGpqanw8vLCzz//XO4PMjY2Fq+++iqKi4vRo0cPNGzYECdOnEBubi769OmDdevWlQukpkyZgri4OAwcONBo3QICAvD8889LbxQzqu32LS0lJQUTJkxAZmYm5syZU+GH/NKlS/Htt9/C2dkZwcHBKCwsRHx8PLRaLV599VW89tprktrCEmyxfXn/imvfwsJCvPDCCzh+/DicnJzQuXNnuLi44Pz588jMzESjRo3w448/okOHDmWuY0v3L2Cbbcx7WPx7xLvvvosdO3bAyckJ3bt3h729PU6fPo379+8jODgY3377LRwcHMqcYzX3sECV0mq1wuDBgwWVSiVs27bNcDw/P1946aWXBJVKJXz88ceiy3v22WcFlUolfPHFF4ZjRUVFwgcffCCoVCph6tSpZfLn5OQI3bp1EwICAoTDhw8bjmdnZwtPPvmkoFKphG+//bbcdXr37i107NhRUKvVprzcGlfb7VvakSNHhNDQUEGlUlXYrvp8KpVK6N+/v3Dz5k3D8b/++ksIDg4W/P39hT///FN0nS3JFttXEHj/im3f1atXCyqVSnjssceE1NRUw3GNRiO88847hrTi4mJDmi3dv4Jgm20sCLyHxbbvzp07BZVKJQwcOFBIS0szHM/KyhKeeOIJQaVSCd98802Zc6zpHmYQVYVffvlFUKlUwqRJk8ql3b17V+jSpYsQGBgo5OTkVFnWiRMnBJVKJTzyyCOCTqcrk1ZYWCgMGDBAUKlUwpUrVwzH9X/A7733XrnykpKSBJVKJfTp00coKioyHL9165agUqmExx9/3IRXWjtqu30FQRDu3LkjzJs3T+jQoYMQEBBgyFfRh/ykSZMElUol/PLLL+XStmzZIqhUKuH111+vsr41wRbbl/ev+PYdOHCgoFKphOPHj5crr6CgQOjZs6egUqmE8+fPG47b0v0rCLbZxryHxbev/n7cs2dPufL27NkjqFQqYdy4cWWOW9M9zInlVYiJiQEAhIWFlUtr3LgxgoODodVqcfjwYdFlDR48GHZ2ZZteoVBg0KBBAIADBw6UO8fYGG+7du2gUqmQmZmJs2fPGo5fuHABABAUFFRlnWpbbbcvAHz++efYunUr2rRpgw0bNiA4OLjCa+Tl5SE+Ph5yuRyDBw8ul/7II49AJpPh4MGD0Ol0VdbZ0mytfQHev2LbV6PRwMfHB+3atUOXLl3Klefg4ABfX18AJXNUANu7fwHba2OA97Ap7xFff/01du3aZUgrrbi4GAAgl8sNx6ztHmYQVYXLly8DAPz9/Y2m+/n5AQAuXbokuiyVSiW6rCtXrph8jv4PuGHDhnj//ffx8MMPIygoCGFhYVi6dCnu3btXZV1rSm23LwC0adMGn376KXbt2oVu3bpVeo3k5GTodDp4e3sbnSPQpEkTNG3aFBqNBteuXauyzpZma+0L8P4V275OTk7YsGED9u7dW26+CFDyYZOcnAwAaN68OQDbu38B22tjgPewKe8RDg4O8Pf3h5OTU5m8ycnJ+PLLLwEAo0ePLnPcmu5hBlFVyMjIAAB4enoaTffw8CiTz5xl5ebmQqPRmHx9/R/w+vXrceDAAahUKnTp0gV///03vv32W4SHh5f51lSbarN99aZNm4Ynn3yywidATLmGqXW2NFtrX4D3rznKAoBVq1ZBo9GgTZs2hknPtnb/ArbXxgDv4eqU9dZbb2HMmDF47LHHkJ6ejrfffhtPPPGE6GuYWufq4hIHVVCr1QBQLkrW0x/X5xNTllKpFFWW/r8ODg7lukYru/7FixcBABMmTMCcOXMM36DS09Mxa9YsnDx5EnPmzMGPP/5YZZ0trTbbV4oHDx6UKcsYR0fHal/HXGytfQHev+Zo319//RU//PAD7OzsMG/ePMhkMgC2d/8CttfGAO9hqe2bl5eHXbt2Gf5tZ2eHGzduIC8vz9DrZG33MIOoKsjlchQXF5f5AzFGELFShH5cV2xZ+sCpqvz/vv6+fftw69YttG/fvsy5np6e+OyzzzB06FAcO3YMiYmJFXbh1pTabF8pxF6jutcxF1trX4D3b3XL2rp1Kz766CMIgoDZs2ejb9++Jl9DzHVqiq21McB7WGpZDg4O+OOPP6BUKnHu3DksXrwYW7ZswcWLF/HTTz/Bzs7O6u5hDudVwcXFBQCQn59vNF0/3FZR5F2dsvT5CwoKDBPsxFzf1dUVKpXK6E3WvHlzBAQEAChZwLO21Wb7SqG/hr4sYwoKCgDA5AXqLMHW2hfg/Su1rOLiYixbtgzvv/8+dDod3nzzTbzwwgtGr2Er9y9ge20M8B6WWpaDgwPc3d3h4uKC3r174/vvv0ezZs3w559/GiatW9s9zCCqCvqx1czMTKPp+jFXfT5zluXq6mrowjTH9fX0EyArutFrUm22rxT6cfiKrmGu65iLrbWvGLx/y5elVqvxyiuv4Ouvv4ZCocCSJUswderUcvls7f4FbK+NxeA9LO7eaty4Mfr37w/gn3lm1nYPM4iqgr6rNSkpyWi6/riYLlkpZemfchB7TmJiIt5++228//77Fdbj9u3bAMo+TVJbart9TeXn5we5XI60tDSjb4B3795FVlYWHB0d0apVK8nXMRdba1/ev6aXdffuXYwbNw4xMTFo1KgRvv/+e4wcOdJoGbZ2/wK218a8h8WXVVBQgCVLlmDGjBmG3qN/088nKyoqAmB99zCDqCroo+Do6OhyadnZ2Th+/DgUCgX69Okjuqz//e9/5YbntFqtYe2MgQMHirp+cnIyLl++jCZNmuA///kPgJJu0l9++QVbt27FjRs3yp2TkpKChIQEKJVK9OrVq8o6W1ptt6+pHB0d0bt3b2i1WkP3cmnR0dEQBAF9+/Y1uqdhTbO19uX9a1r7PnjwAJMnT8bFixfRqlUrbN26FT179qzwGrZ2/wK218a8h8W3r6OjI6KiorBv3z7ExsaWK6+wsBBHjhwB8M+aW9Z2DzOIqkJYWBi8vb0RFxeHTZs2GY5rNBrMnTsXarUa4eHhcHd3N6RptVokJycjOTkZWq3WcLxbt24ICgpCcnIyli9fbpj0ptPpsGDBAty+fRv9+vVDx44dDeeMHj0arq6uiIyMLLNAWU5ODt577z0AwPPPP2+4WVq2bGmY5DhnzhzcvXvXcM6dO3cwY8YM6HQ6TJ48GQ0bNjRnU0lS2+0rxcSJEwEAS5YsKbMOyaVLl/DFF18AKHms3xrYWvvy/jWtfT/99FNcunQJnp6e2Lhxo6hv3rZ0/wK218a8h01r33HjxgEAFi5ciOvXrxuOq9VqzJs3D9euXYOfn1+ZxTit6R7mBsQiHD9+HFOnToVGo0FgYCB8fX1x5swZZGRkICAgABs2bCiz6NfNmzcNK6keOHDAsKItULJ45vjx45GTk4O2bduiffv2+Ouvv3Djxg34+Phgy5Yt5da/iIqKwpw5cyAIArp27YqmTZvixIkTyMnJwUMPPYSvvvqqTMSdnp6O8ePH48aNG2jYsCG6du0KQRAQHx8PjUaDIUOGYPny5aLX7bG02m7ff3vnnXewc+fOSjfI/fDDD/HTTz8ZvhXpdDocP34cWq0WM2bMwMsvv2yGljEPW2tf3r/i2vfq1at47LHHUFxcjMDAQLRt27bCOk6cOBGdO3c2/NuW7l/A9tqY97D49witVovXXnsNsbGxUCgU6N69OxwcHHD+/HncvXsXPj4++P7778sFr9ZyDzOIEuny5ctYtWoV4uPjoVar4evriyFDhmDKlCnlVk2t7AYDgLS0NKxatQqHDx9Gbm4uvLy80L9/f7z00ktlov3STp48iXXr1iEhIQFFRUVo0aIFRo8ejWeffdboSrr379/Ht99+i/379+PmzZtwcHCASqXCk08+iVGjRol6PLQm1Xb7liYmiBIEATt27MCWLVuQlJQER0dH+Pn5YfLkyXj44YcltoLl2Fr78v6tun1/+OEHLFq0SFT9vvjiCzz66KOGf9va/QvYXhvzHhb/HlFcXIytW7ciMjISV65cQXFxMVq2bImwsLAKe+ys5R5mEEVEREQkAedEEREREUnAIIqIiIhIAgZRRERERBIwiCIiIiKSgEEUERERkQQMooiIiIgkYBBFREREJAGDKCIiIiIJGEQRERERSWAdG/cQ1SGXLl3C9u3bcfToUaSnp6OgoABNmjRB+/btMWDAAISHh8PJyanC8+/cuQNXV9dyWy3UpEGDBiEtLQ2ffvopnnzySUllHD16FHv37kVCQgLS0tJQUFAANzc3tGjRAqGhoQgPD4ePj4+Za267jh8/bthYVazBgwdjzZo1FqoREVWFQRSRGa1cuRJfffUViouL4erqipYtW0KhUCAzMxOHDx/G4cOH8e2332L16tUIDAwsc25hYSG++uorrF+/Hrt27arVIKo6UlJS8O677+LMmTMAAIVCAR8fHzRs2BA5OTn4888/kZCQgG+++QbTp0/HK6+8Uss1tj6dOnUyuifmv/n5+dVAbYioIgyiiMwkMjISq1evhrOzMxYtWoSwsDDI5XJDenJyMt577z0kJCRgypQp2Lt3L5o0aWJIz8jIsPlehUuXLmHSpEnIycmBh4cHXn/9dQwdOrRMQJieno7169cjIiICK1euhCAIePXVV2ux1tbniy++KLepKxFZH86JIjKTtWvXAgDmzJmDRx99tEwABQDt2rXDV199haZNmyI7OxsRERG1UU2L0Wg0eO2115CTk4N27dphx44dePLJJ8v1qHl6euLdd9/FvHnzAJS02/Xr12ujykRE1cIgisgM7t27hxs3bgAA/vOf/1SYr0mTJnj44YcBAGfPnq2RutWUiIgI3LhxA/b29vjiiy/QrFmzSvOPGzcOgYGB0Gq12LFjRw3VkojIfDicR2QG9vb//CnFxsYiICCgwryvvfYaJk6ciKZNmxqOTZgwAfHx8YZ/P/LIIwBKApPg4GDD8XPnziEiIgInTpzA33//DWdnZ/j7++OJJ57AqFGjyvV+6R09ehSbN2/GmTNnkJOTA1dXV3Tq1Aljx441XKu6tmzZYqh7+/btRZ0zbdo0XLlypcI6nDhxAhs2bMDp06eRk5ODhg0bokuXLpgwYQJCQkLK5ddPiN+7dy+ysrLw7bff4s8//4RarYavry+GDh2KKVOmwMXFpcx5/v7+AIDvv/8eoaGh5crV/35effVVvPbaa4bjGo0GERERiImJwfXr15GXlwd3d3d069YNzz77LLp37y6qHarrwYMH6Nu3L9RqNVatWoWwsDCj+SZPnowjR45g+vTpeOONNwzH//77b6xfvx4HDx5EWloa7Ozs0LZtWzz22GMYN24cHB0dy5Tz5ZdfYtWqVXjxxRfx/PPPY82aNYiJiUFGRgYaNmyI4OBgvPTSS4Z21XvnnXewc+dOfPTRR3jooYewevVq/PHHH7h79y6aNGmChx56CNOnT69wKNPU+wEA9uzZgx07duDq1avIzMxEgwYNEBAQgCeeeALDhw+HnV3ZvoT09HR88803iI+Px82bNyEIApo3b47Q0FA899xzHGalMtgTRWQGzs7O6NatG4CSD5i3334bJ06cgE6nK5e3WbNm8PPzQ+PGjQ3HVCoVOnXqZPh3YGAgunXrhgYNGhiOffPNNxg7dix27dqF+/fvw9/fH66uroiPj8fcuXPx3HPP4f79++Wu98knn+C5557D/v37odVq0aFDBygUChw+fBivvfYa3njjDWi12mq9/kuXLuHWrVsASp4YE2vIkCF49dVXoVKpyqV99tlnGD9+PKKjo1FYWAiVSgU7OzscOHAAzz33HJYuXVphudu2bcPEiRNx7NgxeHp6okmTJrh69SpWr16NKVOmGP29mKqwsBDPPfccli1bhrNnz6JRo0Zo37498vLysHv3bowbNw7btm2r9nXEcHFxwaOPPgoA2LVrl9E86enpOHbsGABg9OjRhuOnTp3CY489hu+++w43btxAixYt4O3tjQsXLmDJkiUYO3YsMjMzjZZ569YtjBw5Ehs3bgRQMmSdnZ2NvXv34qmnnsKFCxeMnnfx4kU88cQT+OWXX6BUKtGqVSukp6dj+/btePLJJ3H79u1y50i5HxYtWoRZs2YhLi4OMpkM/v7+sLe3R1xcHN566y288847ZfLfuHEDo0aNwoYNG5CamgofHx/4+voiNTUVGzZswBNPPIGLFy8afU1UTwlEZBYXLlwQunTpIqhUKsNPt27dhBdffFFYt26dkJCQIOh0ugrPT01NNZx37dq1Mmn79u0zpK1YsUIoKCgwpB09elQIDQ0VVCqV8NJLL5U577vvvhNUKpUQEBAgbNy4scz19+7da6jvJ598Uua8gQMHCiqVSti6dauo1x4ZGWmo340bN0SdU5ktW7YIKpVK6NGjh/Drr78ajhcXFwt79uwx1Pvf9dPXW6VSCe+8845w7949w3kbN240pP3vf/8rc57++B9//GG0PuPHjxdUKpWwcuVKw7HNmzcLKpVKeOSRR4S0tDTDcY1GI3z88ceCSqUSunfvLmg0GlGv+dixY4Z6pKamijqntBMnTggqlUro1KmTkJubWy79m2++EVQqlfDss88ajt25c0fo1auXoFKphHnz5pU57/r168KTTz5Z7hxBEISVK1ca6jpkyBDh7NmzhrTk5GShX79+gkqlEqZPn17mvLfffttw3tixY4WrV68a0k6fPi107dpVUKlUwqefflrmPCn3Q1JSkqBSqYSgoCDh2LFjZcrbuXOn0KFDB0GlUglnzpwxHH/jjTcElUolvPbaa0JeXp7heGZmpvDUU08JKpVKeP7558u1LdVf7IkiMpOAgABs27atzBBOXl4eDh06hGXLlmHs2LHo27cvPv/8c+Tn55tU9ueffw4AeOqppzBjxowyj7/37t0bq1atAgDExMTg5MmTAICCggJ89dVXAIDXX38d48aNKzN0MXToUHz66acAgM2bN+PmzZsSXnWJv//+2/D/pYcppSgsLMSXX34JAFi4cCFGjBhhSJPJZBg2bBjeeustACW9fkVFReXK6NChAxYuXGjoyZPJZBg3bpxheOnUqVPVqiNQ0vsGAP369YO3t7fhuKOjI9555x307dsXYWFhyMnJMbnswYMHw9/fv8qf0nr06IFWrVqhsLAQ0dHR5cr89ddfAZTthfruu++Qk5ODQYMG4ZNPPkHDhg0NaS1btsSaNWvg6uqKkydP4tChQ0brumzZMgQFBRn+3bZtWzz33HMAgNOnTxs9R6FQYNWqVWjTpo3hWNeuXQ11K32e1PshMTERANCmTZsyQ+IAMHLkSDzzzDMYPnw4CgsLDcf1v9MRI0aUGfJ1d3fH3Llz8dBDD3FZCSqDQRSRGfn5+WHz5s345Zdf8Oqrr6Jr165QKBSG9KysLKxduxYjRozAnTt3RJV57do1pKSkAAAmTZpkNE/Xrl3RtWtXAMCBAwcAACdPnsS9e/dgb2+PcePGGT1v2LBh8PT0hE6nw8GDB8W+zHKKi4urzLNmzZpKAwL9h9+ZM2fw999/w8XFpcKhwREjRsDOzg7p6elGh1cGDBgAmUxW7njbtm0BwOiwp6lat24NANi+fTs2b96Mu3fvGtIcHBzw3XffYdGiRfD09DS57E6dOqFbt25V/vzbqFGjAJQf0vvrr79w+fJlODs7G4b9AOD3338HgDKBSWnu7u7o06cPgJK5fv/m4eFRbr0zoOp27tSpk9EHD4ydJ/V+aNWqFYCSwGjJkiW4du1amXM++OADLFu2DL169TIc05/z2Wef4ffff4dGozGkBQUF4dtvv8W7775rtA5UP3FiOZEFdOzYER07dsRrr72G/Px8nD59GnFxcfj111+RlZWFGzduYMaMGfj555+rLOvq1asAAKVSiXbt2lWYr1OnTjhz5owh4NKf16pVqwoX7pTJZAgICEB6errhPClKz+/KysqCs7NzuTze3t7lPvjz8vJw+fLlMseuXLkCANBqtRUGfwAgl8tRXFyMq1evonPnzmXSPDw8jJ6jXyneHHOinnzySWzfvh1JSUn4+OOPMX/+fHTs2BEhISF46KGH0LNnzzIPHJhC6jpRo0aNwsqVK3HixAncvn0bzZs3B/BPL9SQIUMMPSwPHjxAWloagJIAt6IlN/R59PdTaRUFiPp2NtZLaOp5Uu+HwMBAPP7444iKisL69euxfv16+Pj4ICQkBH379sVDDz1U7u9ixowZOH78OFJSUvDKK6/AwcEBXbt2RZ8+fdC/f3906NChwutT/cQgisjClEol+vTpgz59+mDGjBl47733sGfPHiQkJODChQtGv8mXlpeXBwBVrmBe+sOx9HmlJ6cboy9Xf54UpYeWrly5ghYtWpTLM3LkSIwcObLMMWNbneh7IQoLCyscDirt3r175Y5Vtdq3IAhVllsVV1dX/Pzzz1i/fj12796N69ev4+LFi7h48SK+++47NG3aFG+88QbGjh1b7WuJ5eXlhdDQUMTFxWH37t148cUXodPpsHv3bgBlh/L09weAcoGsMcZ6lUr3sprClPOqcz8sXboUvXv3xrZt2/Dnn38iLS0N27dvx/bt2+Ho6IixY8dizpw5hvulY8eO2LVrF9atW4f//e9/yMnJwfHjx3H8+HEsX74cKpUKH374IXr06GHiK6a6ikEUkRl88MEHOHbsGEaNGoXp06dXmM/JyQnz5883PCmXkpJSZRClD45Kf+gZo//w0OfX/7eqoat/nydFUFAQmjRpgrt372L//v0YNGiQ5LKUSiWAkicUa3r9qIqCK7VabfS4q6srXn/9dbz++uu4fv264QP30KFDyMrKwvvvv49GjRqZbRkJMcaMGYO4uDhERUXhxRdfxJEjR5CZmQlfX1/07NnTkE/fzgAQFRVl9AlJa1Cd+0EmkyE8PBzh4eG4e/cujh8/jvj4eBw6dAhpaWnYsGEDABgWfgWAFi1a4NNPP8X8+fNx/vx5xMfH4+jRozh+/DguX76MF154Ab/99puhl4/qN86JIjKDgoICXL9+3TDHpDKurq6GgKX0ti8V0c8Tyc/PR3JycoX5zp8/D+CfeR368/TrFxlTXFxcbg6JFHK5HM888wwAYO/evYYhGCn0k42vXbtW4XCQIAg4duwYrl27VmZisFT69bUqKisjI6PcsaysLJw8edIwF6pVq1YYO3Ysli1bhkOHDhmWrNAPpdWUhx9+GG5ubkhMTMS1a9cQFRUFoGSor/Q8sYYNG8Ld3R0AkJSUVGF5iYmJ+Ouvv5Cbm2vZildA6v2Ql5eH8+fPG4YhmzRpgqFDh+LDDz/E77//brhf9b8fQRBw8+ZN/PHHHwAAOzs7dO7cGS+88AK+++47REVFwdXVFfn5+di/f79FXzPZDgZRRGagn5h7/vz5Kr8tx8XFIScnB40aNSqzunnpJ+dK94i0adPG8EHy448/Gi3z9OnThhXQ+/XrBwDo3r073NzcUFRUhE2bNhk9b8+ePcjMzIRMJsNDDz1U1cus1AsvvABfX18UFBTgjTfeMKwbVZGioiKjT5H17NkTDRo0wIMHDypsy6ioKEyaNAlDhw4VPUG/Mvo5Xcbm/Zw9e9ZoEDVlyhSMGzcOO3fuLJfm4uKCLl26ADDP/CtTODg4YPjw4QBKAtoDBw5AJpMZJp2XNmDAAADAxo0bjT4ccP/+fUyaNAkjR46s8N6zNKn3w8qVKzFmzBgsWbKkXH47OzvD4pz6309OTg6GDBmC559/HufOnSt3Tps2bQxPYYp5kILqBwZRRGbQp08fDBkyBEDJ0MCCBQvKLRlQUFCAyMhIw0rRM2bMKDOEVnoy9r8DkBkzZgAAfv75Z6xcubJMj8nx48fx+uuvAwAeeughw4rbSqUSU6dOBVDygbJp06Yyb/7R0dH44IMPAABjx44t87i5FM7Ozli3bh2aNWuGpKQkPPHEE/j+++/LLH8AANnZ2fj5558xbNgwQ3DXrl07QxDp7OxsqPeCBQsQGRlZpt6///47PvzwQwAlyzS0bNmyWvUGYFiW4vvvvy/T23fu3DnMmjXL6DlPPPEEAGDVqlX4v//7vzJpJ0+eNPRw9O/fv9r1M5V+7tO3336LvLw8BAcHw8fHp1y+qVOnwtnZGadOncJbb71V5gnDtLQ0TJ06FdnZ2WjQoEGlk7otSer9MGLECMhkMhw8eBDffPNNmQVlb926ZdjrUv/7ady4seGLxHvvvVfmPiguLsamTZtw+fJls3zhoLqDc6KIzOSzzz6Ds7MzfvnlF0RERCAiIgLe3t5o2rQpCgoKDEMNCoUCb775Jp599tky5zdq1Ag+Pj5IS0vDK6+8grZt22LGjBno168fhg4dihs3buDzzz/H6tWr8eOPP6JNmza4e/eu4empXr16YenSpWWGbKZMmYKbN29iy5YtmD9/Pr788ku0aNECd+7cMfSuDBkyBHPnzjVLG/j5+WHbtm346KOPcPDgQSxevBj//e9/4eXlhaZNmyI3NxdpaWmGb//NmjXD5MmTMXHixDI9cS+++CJSU1OxdetWvPfee1i6dCl8fX2Rnp5uqHe3bt0M61xV1/Tp03H48GFkZmbi8ccfh5+fn+F31qJFC4wZMwaRkZFlzpk4cSKOHDmC//u//8OLL74IDw8PeHh4IDs72/A7GTRoEJ588kmT6/PvtcAqs3LlynLLBXTq1AkqlcowYbz0hPLSWrVqhRUrVmDmzJnYvXs3oqOj4efnB61Waxg+c3Z2xtdff13t9b+qQ8r90KlTJ7zxxhv4/PPP8dlnn+Hrr7+Gr68v8vPzkZqaiqKiIrRs2bLMquXz58/HU089hcuXL2P48OHw9fVFgwYNcOvWLWRnZwMAZs2axbWiyIBBFJGZODg4YPHixRg3bhz27t2L48ePIz09HZcuXYJSqUSbNm3Qt29fhIeHG+Yr/dsXX3yBBQsW4K+//sK1a9cMmxoDJfvMhYSE4Mcff8TJkydx6dIlNGzYECEhIRg5cqRhrZzSZDIZPvroIzz88MPYsmULEhIS8Ndff6Fx48YYOHAgwsPDDRsim0vz5s2xbt06XLhwAVFRUThx4gRu3LiB9PR0uLq6ok2bNvjPf/5jWIzS2JNaMpkMn3zyCYYMGYKffvrJUG9HR0d06dIFw4cPx1NPPSU60KhKx44dsX37dnz11Vc4duwYrl69Ci8vLzz//PN4+eWX8e2335Y7Ry6XY/Xq1diyZQt+++03JCcnG34nffv2xYgRIwy9IabSz28To6CgwOjxMWPGYNGiRXBxcal0Ynv//v2xZ88e/PDDDzh8+DBSUlKg0+ng4+ODPn364Pnnnzf6tGVNkno/vPTSS/Dz88PWrVtx4cIFXL58GU5OTujYsSPCwsIwYcKEMj3AHh4e2L59O7777jscPnwYqampuH37Npo2bYrHHnsM48ePN7o+F9VfMsEcz/oSERER1TOcE0VEREQkAYMoIiIiIgkYRBERERFJwCCKiIiISAIGUUREREQSMIgiIiIikoBBFBEREZEEDKKIiIiIJGAQRURERCQBgygiIiIiCRhEEREREUnAIIqIiIhIAgZRRERERBL8P6HuYaBWQ+gSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(x='genus_evenness',\n",
    "                y='median_mmNorm_DTAPHib',\n",
    "                data=stool_meta.query(\"VisitCode == 'V5'\"))\n",
    "plt.xlabel('Stool Genus Evenness')\n",
    "plt.ylabel('Meidan DTAPHib Titer')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387b8bb0-f5f2-469c-8fd7-f4b326bbb393",
   "metadata": {},
   "source": [
    "## Combined analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2fe9cb64-807e-4fa6-b0e1-3770cc2ce54f",
   "metadata": {},
   "outputs": [],
   "source": [
    "titer_stats = pd.concat([nasal_titer_stats, stool_titer_stats]).sort_values('p_value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "0052e66a-f697-4c53-b7f7-ee74b079615d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.454801</td>\n",
       "      <td>0.000906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.387002</td>\n",
       "      <td>0.003841</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>341</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.556600</td>\n",
       "      <td>0.010806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>303</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.478974</td>\n",
       "      <td>0.013302</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>339</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>ko_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.522001</td>\n",
       "      <td>0.018233</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>328</th>\n",
       "      <td>stool</td>\n",
       "      <td>V2</td>\n",
       "      <td>ko_MDS2</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.459145</td>\n",
       "      <td>0.018299</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>144</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V8</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.329268</td>\n",
       "      <td>0.019548</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.313894</td>\n",
       "      <td>0.02081</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>150</th>\n",
       "      <td>stool</td>\n",
       "      <td>V4</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>14.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.604396</td>\n",
       "      <td>0.022057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.352104</td>\n",
       "      <td>0.032579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>349</th>\n",
       "      <td>stool</td>\n",
       "      <td>V3</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>20.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.470677</td>\n",
       "      <td>0.036215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.342579</td>\n",
       "      <td>0.037933</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>genus_evenness</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.450649</td>\n",
       "      <td>0.040345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>97</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>genus_shannon</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.448052</td>\n",
       "      <td>0.041656</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V6</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>52.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.293825</td>\n",
       "      <td>0.045005</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>stool</td>\n",
       "      <td>V5</td>\n",
       "      <td>ko_evenness</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>32.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.354106</td>\n",
       "      <td>0.046769</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>110</th>\n",
       "      <td>stool</td>\n",
       "      <td>V1</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>21.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.433766</td>\n",
       "      <td>0.049467</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint      div_metric          titer_measure N_LVR N_NVR  \\\n",
       "145     nasal        V8            MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "95      nasal        V9   simpson_e_div      median_mmNorm_PCV   NaN   NaN   \n",
       "341     stool        V3     ko_richness      median_mmNorm_PCV   NaN   NaN   \n",
       "303     stool        V2  genus_evenness          median_mmNorm   NaN   NaN   \n",
       "339     stool        V3     ko_richness          median_mmNorm   NaN   NaN   \n",
       "328     stool        V2         ko_MDS2  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "144     nasal        V8            MDS1          median_mmNorm   NaN   NaN   \n",
       "93      nasal        V9   simpson_e_div          median_mmNorm   NaN   NaN   \n",
       "150     stool        V4  genus_richness          median_mmNorm   NaN   NaN   \n",
       "180     stool        V9  genus_richness          median_mmNorm   NaN   NaN   \n",
       "349     stool        V3     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "199     stool        V9     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "95      stool        V1  genus_evenness      median_mmNorm_PCV   NaN   NaN   \n",
       "97      stool        V1   genus_shannon  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "129     nasal        V6            MDS1          median_mmNorm   NaN   NaN   \n",
       "43      stool        V5     ko_evenness  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "110     stool        V1     kraken_MDS1      median_mmNorm_PCV   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value  \n",
       "145  50.0  NaN  0.454801  0.000906  \n",
       "95   55.0  NaN -0.387002  0.003841  \n",
       "341  20.0  NaN -0.556600  0.010806  \n",
       "303  26.0  NaN -0.478974  0.013302  \n",
       "339  20.0  NaN -0.522001  0.018233  \n",
       "328  26.0  NaN -0.459145  0.018299  \n",
       "144  50.0  NaN  0.329268  0.019548  \n",
       "93   55.0  NaN -0.313894   0.02081  \n",
       "150  14.0  NaN  0.604396  0.022057  \n",
       "180  37.0  NaN  0.352104  0.032579  \n",
       "349  20.0  NaN  0.470677  0.036215  \n",
       "199  37.0  NaN  0.342579  0.037933  \n",
       "95   21.0  NaN -0.450649  0.040345  \n",
       "97   21.0  NaN  0.448052  0.041656  \n",
       "129  52.0  NaN  0.293825  0.045005  \n",
       "43   32.0  NaN -0.354106  0.046769  \n",
       "110  21.0  NaN  0.433766  0.049467  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "titer_stats.query('p_value < .05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "92462504-4b5b-4d9d-9df8-0341a760c2c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>95</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm_PCV</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.387002</td>\n",
       "      <td>0.003841</td>\n",
       "      <td>0.691446</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.342579</td>\n",
       "      <td>0.037933</td>\n",
       "      <td>0.976296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>93</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V9</td>\n",
       "      <td>simpson_e_div</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>55.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.313894</td>\n",
       "      <td>0.02081</td>\n",
       "      <td>0.976296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.352104</td>\n",
       "      <td>0.032579</td>\n",
       "      <td>0.976296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>stool</td>\n",
       "      <td>V5</td>\n",
       "      <td>ko_evenness</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>32.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.354106</td>\n",
       "      <td>0.046769</td>\n",
       "      <td>0.976296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>129</th>\n",
       "      <td>nasal</td>\n",
       "      <td>V6</td>\n",
       "      <td>MDS1</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>52.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.293825</td>\n",
       "      <td>0.045005</td>\n",
       "      <td>0.976296</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint      div_metric          titer_measure N_LVR N_NVR  \\\n",
       "95      nasal        V9   simpson_e_div      median_mmNorm_PCV   NaN   NaN   \n",
       "199     stool        V9     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "93      nasal        V9   simpson_e_div          median_mmNorm   NaN   NaN   \n",
       "180     stool        V9  genus_richness          median_mmNorm   NaN   NaN   \n",
       "43      stool        V5     ko_evenness  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "129     nasal        V6            MDS1          median_mmNorm   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value     p_adj  \n",
       "95   55.0  NaN -0.387002  0.003841  0.691446  \n",
       "199  37.0  NaN  0.342579  0.037933  0.976296  \n",
       "93   55.0  NaN -0.313894   0.02081  0.976296  \n",
       "180  37.0  NaN  0.352104  0.032579  0.976296  \n",
       "43   32.0  NaN -0.354106  0.046769  0.976296  \n",
       "129  52.0  NaN  0.293825  0.045005  0.976296  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vaccine_times = ['V5', 'V6', 'V7', 'V9']\n",
    "titer_stats_vax = titer_stats.query('timepoint in @vaccine_times').sort_values(['timepoint', 'div_metric', 'titer_measure'])\n",
    "titer_stats_vax['p_adj'] = p_adjust(titer_stats_vax['p_value'])\n",
    "titer_stats_vax = titer_stats_vax.sort_values('p_adj')\n",
    "titer_stats_vax.query('p_value < .05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "da67fd63-32e4-47ec-ad7d-9925934c4281",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>body_site</th>\n",
       "      <th>timepoint</th>\n",
       "      <th>div_metric</th>\n",
       "      <th>titer_measure</th>\n",
       "      <th>N_LVR</th>\n",
       "      <th>N_NVR</th>\n",
       "      <th>N</th>\n",
       "      <th>stat</th>\n",
       "      <th>R</th>\n",
       "      <th>p_value</th>\n",
       "      <th>p_adj</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>kraken_MDS1</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.342579</td>\n",
       "      <td>0.037933</td>\n",
       "      <td>0.942978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>stool</td>\n",
       "      <td>V9</td>\n",
       "      <td>genus_richness</td>\n",
       "      <td>median_mmNorm</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>37.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.352104</td>\n",
       "      <td>0.032579</td>\n",
       "      <td>0.942978</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>stool</td>\n",
       "      <td>V5</td>\n",
       "      <td>ko_evenness</td>\n",
       "      <td>median_mmNorm_DTAPHib</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>32.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.354106</td>\n",
       "      <td>0.046769</td>\n",
       "      <td>0.942978</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    body_site timepoint      div_metric          titer_measure N_LVR N_NVR  \\\n",
       "199     stool        V9     kraken_MDS1  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "180     stool        V9  genus_richness          median_mmNorm   NaN   NaN   \n",
       "43      stool        V5     ko_evenness  median_mmNorm_DTAPHib   NaN   NaN   \n",
       "\n",
       "        N stat         R   p_value     p_adj  \n",
       "199  37.0  NaN  0.342579  0.037933  0.942978  \n",
       "180  37.0  NaN  0.352104  0.032579  0.942978  \n",
       "43   32.0  NaN -0.354106  0.046769  0.942978  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vaccine_times = ['V5', 'V6', 'V7', 'V9']\n",
    "titer_stats_vax_stool = titer_stats.query('timepoint in @vaccine_times and body_site == \"stool\"').sort_values(['timepoint', 'div_metric', 'titer_measure'])\n",
    "titer_stats_vax_stool['p_adj'] = p_adjust(titer_stats_vax_stool['p_value'])\n",
    "titer_stats_vax_stool = titer_stats_vax_stool.sort_values('p_adj')\n",
    "titer_stats_vax_stool.query('p_value < .05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a58bd5e4-4a64-474f-b89e-a194f2b2083d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
